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CHAPTER  1 
INTRODUCTION 

1.1   IDEA  OF  DECISION  MAKING  IN  MANAGEMENT 

The  administration  of  a  modern  business  enterprise  has  become  an 
enormously  complex  undertaking.   During  the  past  few  years  there  has  been 
an  increasing  tendency  to  turn  to  quantitative  techniques  and  models  as  a 
potential  means  for  solving  the  problems  that  arise  in  such  an  enterprise. 

Engineering  has  been  defined  as  concerned  with  the  design,  improve- 
ment and  installation  of  integrated  systems  of  men,  machines  and  materials 
for  the  service  of  society.   Every  working  day,  the  typical  executive  of 
a  modern  industrial  organization  makes  a  number  of  complex  decisions  in 
order  to  optimize  his  company's  performance.   This  emphasizes  the  impor- 
tance of  quantitative  techniques  as  a  useful  means  for  decision  making  in 
industrial  management  systems. 

In  any  problem  solving  situation,  there  are  variables  or  factors 
which  influence  the  outcome  of  whatever  decision  is  made.  These  variables 
can  be  classified  as  those  which  the  decision  maker  controls,  called  the 
control  variables,  and  a  class  of  those  which  he  cannot  control,  called 
the  state  variables.  After  identifying  the  control  and  state  variables, 
they  should  be  combined  in  some  logical  manner  so  that  they  form  a  model 
of  the  problem.   The  object  of  the  decision  maker  is  not  to  construct  a 
model  as  close  as  possible  to  the  reality  of  the.  problem,  but  rather  a 
simplest  model  that  predicts  the  outcomes  reasonably  well.  Next  step  is 
to  develop  a  measure  of  effectiveness,  called  the  objective  function,  to 
predict  the  behavior  of  the  model.  A  model  is  then  solved  for  different 


values  of  the  control  variables.  These  will  be  called  feasible  solutions. 

In  general,  decision  making  can  be  described  as  a  process  whereby 
management  when  confronted  with  a  problem,  selects  a  specific  course  of 
action,  called  the  optimal  policy,  from  a  set  of  feasible  solutions. 

Many  of  the  mathematical  models  in  engineering,  physical  sciences  and 
other  disciplines  involve  non-linear  differential  equations  of  the  two 
point  boundary  value  type.  Unfortunately,  no  general  analytical  method 
exists  for  solving  them.   Several  kinds  of  non-linear  differential 
equations  have  been  solved  analytically,  but  the  solution  of  each  has 
required  a  method  unique  to  that  type.   Various  methods  have  been  used 
to  solve  non-linear  differential  equations  numerically.  Among  them  are 
graphical  methods,  methods  based  upon  successive  approximations  and  methods 
based  upon  iterative  procedures. 

1.2   PURPOSE  OF  THIS  STUDY 

Industrial  Engineers  work  with  a  wide  variety  of  optimization  problems. 
For  this  reason  they  should  be  familiar  with  the  most  efficient  techniques 
for  solving  the  decision-making  problems.   Because  of  the  relatively  recent 
origin  of  operations  research,  more  efficient  techniques  are  not  needed  in 
most  cases. 

The  purpose  of  this  research  is  to  study  the  effectiveness  of  a  recently 
developed  method,  quasilinearization,  in  solving  industrial  management 
problems  which  involve  non-linear  differential  equations. 

More  specifically,  the  object  of  this  work  is  to  investigate  the  compu- 
tational features  of  this  technique  with  respect  to  different  problems. 
The  second  object  is  to  provide  the  systems  analysts  a  new  tool  for 


optimization. 

Other  computational  techniques,  such  as  the  gradient  technique,  the 
second  variation  method,  and  invariant  imbedding  can  also  be  used  for 
solving  the  problems  with  non-linear  differential  equations,  but  consid- 
ering the  object  of  this  study,  they  will  not  be  discussed  here. 


CHAPTER  2 
SOLUTION  OF  TWO  POINT  BOUNDARY  VALUE  PROBLEMS 

2.1   INTRODUCTION 

The  mathematical  formulation  of  many  problems  in  science  and  engin- 
eering leads  to  differential  equations.   Problems  in  which  the  conditions 
to  be  satisfied  by  the  solution  of  a  differential  equation  of  order  two 
or  greater  may  be  specified  at  both  ends  of  an  interval  are  known  as  two 
point  boundary  value  problems.   If  the  conditions  are  specified  at  more 
than  two  points  in  the  interval,  the  problems  are  known  as  multi-point 
boundary  value  problems.   The  latter  type  of  problems  do  not  appear  very 
often  in  engineering  models. 

Initial  value  problems  are  those  in  which  all  conditions  are  imposed 
at  one  point.  This  may  be  the  initial  or  final  point  of  that  interval. 

Consider  a  system  of  differential  equations 


dx 


(i) 


ar =  8(x,y) 


If  the  conditions  for  both  x  and  y  aie  given  at  the  same  point v 


xl(tl)  =  Xl         X2(tl)  =  X2 


X;L(tf)  =  x^         x2(tf)  -  x2  (2) 


the  problem  is  called  the  initial  value  problem.   However,  if  the  con- 
ditions for  both  x  and  y  are  not  at  the  same  point, 


X1(t±)     =    HLl 

/       ^             1 

x2(tf)   ■  x2 

or 

X    (tf)    =    JE. 

x2(t.)   =  x2 

(3) 


the  problem  is  called  two  point  boundary  value  problem. 

A  higher  order  differential  equation  can  always  be  replaced  by  a  set 
of  first  order  differential  equations  by  introducing  auxiliary  variables  [9]. 
For  this  reason,  only  first  order  differential  equations  will  be  discussed 
throughout  this  work. 

2.2   NUMERICAL  SOLUTION  OF  INITIAL  VALUE  PROBLEMS 

Since  all  of  this  work  will  be  based  on  numerical  methods  of  obtaining 
solutions  of  ordinary  differential  equations,  a  best  known  and  most  fre- 
quently used  scheme  for  solving  initial  value  problems,  Runge-Kutta  method, 
is  discussed  here.  ' 

In  this  method,  the  increments  of  the  functions  are  calculated  once 


for  all  by  means  of  a  definite  set  of  formulas ,  and  the  calculations  for 
the  first  increment  are  exactly  same  as  for  any  other  increment.   These 
processes  are  self-starting.   Advantage  of  this  method  is  the  independent 
choice  of  the  size  of  the  step,  which  may  be  increased  to  speed  up  the 
progression  or  decreased  to  lower  truncation  errors  without  recalculation 
of  previous  data. 

The  fourth-order  formulas  for  the  Runge-Kutta  method,  [9],  for 
Equations  (1)  are 


X1(W  "  Xl(tk5  +  I  S  +  2m2  +  2m3  +  V 
"2(W  =x2(tk)+i(nl+2n2  +  2n3  +  V 


where 


n^fCx^),  x2(tk),  tk)  At 

n2  =  f (xl(tk)  +  -|,  VV+1'  tk  +  Al)  At 


n3=f(x1(tk)+-f,X2(tk)+-f,  t^^at 


i4  •  f (x1(tk)  +  m3>  *2(tk)  +  n3>  \  *   At)  At 


n2  -  g(xx(tk)  +  -i,  x2(tk)  +  -|,  tfc  +  -§)  At 
n3  =  g(xl(tk)  +  \   x2(tk)  *\t±*  %  At 


n4  =  g(x1(tv.)  +  m3,  x2(tk)  +  n3,  tfc  +  At)  At. 


w 


(5) 


READ  AT,  tf,  x(t  )  and 
i&J , 


D 


Compute  m  ,  m  ,  m  ,  m. 


Fig.  1.   Computer  Logic  Diagram  for  Eunge-Kutta  Method 


Knowing  the  initial  values  of  x^t  )  ,  x,(t.),  and  step  size  At, 
values  of  x  (t.  +  At)  and  x  (t .  +  At)  can  be  calculated  using  the  above 
formulas.   Similarly  ^(t.  +  2At)  and  x2(t±  +  2At)  can  be  calculated 
using  x  (t.  +  At)  and  x  (t.  +  At).   Hence  incrementing  t  everytime  by 
At,  the  final  values,  x  (t  )  and  x2(tf)  can  be  calculated. 

The  truncation  error  in  this  method  is  0(At  ).  A  simplified  compu- 
tational scheme  is  shown  in  Figure  1. 


2.3   DIFFICULTIES  IN  TWO  POINT  BOUNDARY  VALUE  PROBLEMS 

The  numerical  solution  of  any  ordinary  differential  equation  requires 
the  knowledge  of  initial  values  of  all  the  variables.  Starting  with  the 
initial  values,  the  solutions  are  constructed  step  by  step  in  small 
intervals  of  the  variables.   Because  of  this  nature,  they  are  also  called 
the  marching  techniques. 

In  an  initial  value  problem,  all  the  initial  (or  final)  values  are 
known.  Hence  the  solution  is  relatively  easy.   In  a  two  point  boundary 
value  problem,  some  of  the  initial  (or  final)  values  are  unknown.   Hence, 
the  numerical  techniques,  like  Runge-Kutta  method,  cannot  be  applied 
directly.   For  this  reason,  this  type  of  problems  are  very  difficult  to 
solve. 

In  general,  the  procedure  for  solving  this  type  of  problems  is  to 
assume  the  missing  initial  (or  final)  conditions  and  solve  for  all  the 
grid  points  and  then  compare  the  values  of  the  calculated  and  given  final 
(or  initial)  conditions.   If  they  are  not  the  same  within  allowable  error, 
a  new  set  of  missing  initial  (or  final)  values  is  assumed  and  the  same 


procedure  is  repeated.   By  this  trial  and  error  procedure,  a  suitable  set 
of  initial  (or  final)  conditions  can  be  determined. 

This  procedure  becomes  very  tedious  if  the  problem  has  many  differen- 
tial equations  and  is  very  complex  in  nature.  The  relatively  slow  conver- 
gence during  the  process  of  numerical  solution  can  make  the  generally  used 
trial  and  error  procedure  impractical. 

Unfortunately,  most  of  the  mathematical  models  in  quantitative  analysis 
are  very  complex  having  many  variables.   Problems  of  this  type  are  most 
subtle  and  difficult  and  are  not  well  suited  for  modern  digital  computers. 
There  is  no  general  proof  of  existence  and  uniqueness  of  solutions  to 
problems  of  this  type. 

2.4   SUPERPOSITION  PRINCIPLE 

A  two  point  boundary  value  problem  is  not  too  difficult  if  the  per- 
formance equations  are  linear.   This  is  because  of  the  fact  that  super- 
position principle  is  applicable  to  linear  differential  equations. 

Consider  the  following  two  simultaneous  first  order  linear  differential 
equations 


dt 


=  a^t)  +  b1(t)  X;L  +  cx(t)  x2  (6) 


dx 

-7—  -  a  (t)  +  b  (t)  x.  +  c  (t)  x„  (7) 


*1^!  =  xl   and   X2^P  =  x2  ^ 


10 


t.   X.. ,  x  are  known  constants  at  the  initial  and  final  values  of  t 
respectively. 


x_  (t.)  =  1  and  x-  (t.)  =  0,  Equations  (6)  and  (7)  can  be  solved  numerically 

to  obtain  a  set  of  particular  solutions,  x1  (t)  and  x„  (t) ,  t,  <  t  <  t_. 

Ip         2p      i  -   —  t 

Two  sets  of  non-trivial  homogeneous  solutions,  x..  .,  (t) ,  x_  1u(t)  and 

l,±n      2,  in 

x.  _,  (t) ,  x    (t),  can  be  obtained  with  any  two  different  sets  of  arbitrarily 

assumed  initial  conditions,  say  x,  n_(t.)  ■  1,  x.  .,  (t.)  =  0  and  x.  ^,(t.)  =  0, 

l,ln  i       2,  In  a.  l,2n  i 

x«  _,  (t.)  =  1,  from  the  homogeneous  equations  of  Equations  (6)  and  (7). 

The  homogeneous  equations  are  obtained  by  setting  the  constant  terms 
equal  to  zero. 


dx 

3g±  =  b1(t)  Xj  +  c^t)  x2  (9) 


dx 

-£*■  =  b2(t)  Xj  +  c2(t)  x2  (10) 


It  is  important  to  note  that  the  particular  and  homogeneous  solutions 
can  be  obtained  numerically  using  a  step  by  step  integration  method,  like 
the  Runge-Kutta  method.   The  reader  is  referred  to  Ince  [7]  and  Lee  [9] 
for  detailed  discussion. 

The  superposition  principle  states  that  because  of  the  additive 
property  of  the  solution  of  a  linear  system,  the  general  solution  of 
Equations  (6)  and  (7)  is, 


xx(t)  =  x   (t)  +  A^  lh(t)  +  A2x1  2h(t)  (11) 


x2(t)  =  x2  (t)  +  Axx2  lh(t)  +  A2x2  2h(t)  (12) 


where  A.  and  A.  are  integration  constants. 

A.,    and  A_  can  be  obtained  by  substitv 
Equation  (8),  into  Equations  (11)  and  (12)  with  the  results  of  the  partic- 
ular and  homogeneous  solutions.   Once  the  values  of  A_  and  A.  are  known, 
the  right  hand  sides  of  Equations  (11)  and  (12)  are  completely  known. 
This  gives  the  solution  for  x. (t)  and  x-(t)  at  all  the  grid  points. 

This  approach  can  be  generalized  to  a  set  of  n  simultaneous  first 
order  linear  differential  equations 


dx 

j+  =  g..(Xl,  x2,  .  .  .,  xn,  t)        i  -  1,  2,  ...  a    (13) 


x (tf)  =  x*  J  -  X,   2,  ...  a    (14) 


x^(t  )  =  x]<^  k  =  m  +  1,  m  +  2,  .  .  .  n    (15) 


The  general  solution  by  the  superposition  principle  is 


*i(t)  =  xip(t)  +  Z  Vi  kn(t)        1  -  1,  2,  ...  a    (16) 

k=l     ' 


In  this  general  case,  we  have  to  assume  n  initial  conditions, 

x.  (t.)  "  x.  ,  for  the  particular  solution,  and  n  sets  of  initial  conditions 
ip  i     ip 


x    (t  )  =  x.  ,  ,  for  n  sets  of  homogeneous  solutions.   Integration 
constants  A,  are  determined  from  the  n  known  boundary  conditions  and  the 
assumed  and  computed  boundary  conditions  for  the  particular  and  homo- 
geneous solutions . 

Usually  n  sets  of  homogeneous  solutions  are  required  to  obtain  the 
general  solution.   However,  if  the  assumed  initial  values  for  the  partic- 
ular solution  are  properly  selected,  only  m  sets  of  homogeneous  solutions 
need  be  obtained. 

Consider  the  Equations  (6)  and  (7),  their  general  solution  is  given 
by  Equations  (11)  and  (12).   Suppose  the  initial  values  for  the  particular 
solution  are  chosen  as 


xlp(t.)  =  x°  x2p(ti)  =  0 


and  the  initial  values  for  the  homogeneous  solutions  are  given  as  before, 
then  Equations  (11)  and  (12)  at  the  initial  time  t  ,  reduce  to 


«l(ti>  =  VV  +  Vl.ll/V  +  A2xl,2h(ti) 


x°  =  x°  +  Aj       1  +  A2       0  (17) 


A.   ■  x     -  x1   =  0. 


This  shows  how  to  select  the  appropriate  initial  conditions  so  as  to 
reduce  the  set  of  homogeneous  solutions  needed  from  n  to  m.   For  further 


13 


discussion,  the  reader  is  referred  to  Lee  [9]. 

The  procedure  can  be  divided  into  essentially  two  steps.   First,  the 
problem  is  converted  into  initial  value  problems  and  these  problems  are 
solved  numerically.   Then,  the  integration  constants  are  obtained  by 
solving  a  set  of  algebraic  equations.   Combination  of  these  results  is  the 
general  solution  of  the  original  problem. 


14 


CHAPTER  3 
QUASILINEARIZATION 

3.1   INTRODUCTION 

The  advantages  of  superposition  principle  in  solving  two  point  bound- 
ary value  problem  lead  to  the  idea  of  linearizing  the  non- linear  differen- 
tial equations  so  that  the  superposition  principle  can  be  applied.  This 
is  the  basic  concept  of  quasilinearization. 

Quasilinearization  technique  was  developed  by  Bellman  [  1 ]  and 
Kalaba  [  8 ]  and  applied  extensively  to  chemical  engineering  problems  by 
Lee  [q  ,  10,11]  in  obtaining  numerical  solutions  of  certain  classes  of 
non-linear  ordinary  differential  equations  of  the  boundary  value  type 
encountered  in  chemical  engineering,  optimization,  the  boundary  layer 
theory  and  in  control  problems. 

This  technique  essentially  linearizes  the  set  of  non-linear  differ- 
ential equations.   Conceptually,  this  method  is  very  close  to  Newton 
Raphson  method  of  finding  roots  of  an  equation;  however,  since  the 
unknowns  to  be  determined  in  this  method  are  functions  and  not  fixed 
valued  roots  as  in  Newton  Raphson  method,  both  the  computational  and 
theoritical  aspects  are  much  more  complicated. 

In  addition  to  linearizing  the  non-linear  equations,  the  quasilinear- 
ization technique  provides  a  sequence  of  functions  which  in  general  con- 
verges rather  rapidly  to  the  solution  of  the  original  non-linear  equations. 
Usually,  the  latter  is  more  important.  A  rough  initial  approximation  for 
the  unknown  function  can  lead  to  the  solution  of  the  origina]  equation 
through  a  sequence  of  functions.   In  general,  for  most  practical  problems, 
this  rough  initial  approximation  can  be  obtained  from  engineering 
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experiences  and  intuitions. 

3.2   COMPUTATIONAL  PROCEDURE 

In  many  operations  research  techniques,  the  verbal  description  of  the 
algorithm  is  far  more  difficult  than  the  algorithm  itself.   Hence  the 
logic  will  be  developed  and  explained  with  an  illustration. 

Consider  a  set  of  nonlinear  differential  equations 


fc-ffe.y) 


(18) 


dt 


=  g(x,y) 


with  boundary  values 


x(t.)  =  x°       and       y(tf)  =  y1  (19) 


Using  the  Taylor  series  expansion  f (x,y)  and  g(x,y)  can  be  linearized 
around  x  =  a  and  y  ■  b  as  follows: 


f (x,y)  =  f(a,b)  +  (x  -  a)  fa(a,t)  +  (y  -  b)  ^(b.t) 


g(x,y)  =  g(a,b)  +  (x  -  a)  ga(a,t)  +  (y  -  b)  gb(b,t)  (20) 


which  is  the  Taylor  series  with  second  and  higher  terms  omitted.   Symbol 
f  (a,t)  represents  the  partial  derivative  of  f  with  respect  to  x  at  x  ■  a. 
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From  Equations  (18)  and  (20) ,  we  obtain 


j£  =  f (a,b)  +  (x  -  a)  fa(a,t)  +  (y  -  b)  ^(b.t) 


|2-  =  g(a,b)  +  (x  -  a)  8a(a,t)  +  (y  -  b)  gb(b,t)  (21) 


Since  a  and  b  are  known  functions  of  t,  Equations  (21)  are  linear 
differential  equations  with  variable  coefficients.   The  boundary  conditions 
for  Equations  (21)  are  given  by  Equations  (19) . 

A  recurrence  relation  can  now  be  established.   Choose  an  initial 
approximation  for  a  and  b,  say  a  =  x-  and  b  =  y_.   Substituting  these 
approximations  into  Equations  (21) ,  it  is  possible  to  solve  these  first 
order  linear  differential  equations  for  x  and  y  using  a  step  by  step 
integration  method  and  the  boundary  conditions  given  by  Equation  (19). 
Call  this  new  solution  of  x  and  y  as  x..  and  y-  .  Now  using  x..  and  y_ ,  it 
is  possible  to  find  improved  values  of  x  and  y.   Call  these  improved 
functions  x-  and  y  .   Next  using  x  and  y„,  x,  and  y  can  be  determined. 
This  iterative  procedure  is  continued  until  the  desired  accuracy  is 
obtained. 

The  recurrence  relation  can  be  written  as 

dx 

—  =  f  (xo>  yQ)  +  Uj  -  xq)  fx  (xo,  yQ)  +  (7l  -  yQ)  f   (x   y  )     (22) 
o  'o 

*1 

dT"  g(xo'  V  +  (X1  '  Xo}  8x  (xo'  V  +  (yl  -  yo}  gy  (V  yo) 
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'2  =  f(xr   yi)  +   (x2  -  «x)    f^,   yx)   +    (y2  -  yx)    f^,   y^  (23) 


dt 


J^  =  g(Xl,   7X)   +    U2  "  xx)   gXi(xr   7l)  +    (y2  -  yx)    8yi(xr   y^ 


dxN+l 

-dt~"  =   f  (V   V  +   (XN+1  "  XN)    fxH(xN'   V  +    (yN+l  "  VfyN(W    (2A) 


dy 

"IT"  =  8(V   V  +    (XN+1  "  V    «^V   V   +   (yN+l  "  V«JB  W 


The  boundary  conditions  given  in  Equation  (19)  are  used  in  solving  Equations 
(22)  through  (24). 

For  a  number  of  problems  Equations  (22)  through  (24)  have  been  proved 
to  converge  monotomically  to  the  solution  of  Equation  (18) .  The  convergence 
rate  is  quadratic  in  the  sense  that  each  iteration  approximately  doubles 
the  number  of  digits  of  accuracy. 

3 . 3   SUMMARY 

The  procedure  can  be  summarized  in  the  following  steps. 

1.  The  nth  order  non-linear  ordinary  differential  equations  are  first 
converted  into  a  system  of  simultaneous  first  order  ordinary  differ- 
ential equations. 
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2.  This  set  of  equations  is  then  linearized  using  Equation  (21). 

3.  The  recurrence  relation  for  the  set  of  linearized  first  order  differ- 
ential equations  is  constructed  using  Equations  (22)  through  (24) . 

4.  The  appropriate  initial  approximation,  x   (t) ,  is  assumed  for  each 
unknown  dependent  variable  as  a  function  of  independent  variable  t. 

5.  The  results  of  the  first  iteration,  x.  ..(t),  can  be  obtained  by  substi- 
tuting x  (t)  into  the  recurrence  relation  and  using  the  superposition 
principle. 

6.  A  further  improved  solution  is  obtained  by  repeating  Step  5.   The 
procedure  is  continued  until  the  solution  converges  to  the  desired 
accuracy. 

3.4   DISCUSSION 

The  main  advantage  of  this  technique  is  that  if  the  procedure  conver- 
ges, it  converges  quadratically  to  the  solution  of  the  original  equation. 
Quadratic  convergence  means  that  the  error  in  the  (n+l)st  iteration  tends 
to  be  proportional  to  the  square  of  the  error  in  the  nth  iteration.   All 
the  computational  features  of  Newton-Raphson  technique  are  retained  in 
this  technique. 

In  spite  of  all  the  advantages,  this  technique  also  has  its  difficulties. 
There  are  two  main  difficulties.  The  first  difficulty  arises  from  the  fact 
that  in  using  the  superposition  principle,  a  set  of  algebraic  equations  must 
be  solved.   Thus  the  ill-conditioning  phenomenon  in  solving  a  set  of  linear 
algebraic  equations  can  make  the  superposition  principle  useless.  Another 
difficulty  is  the  convergence  problem.   If  the  initial  approximation  is  not 


within  the  interval  of  convergence  a  solution  cannot  be  obtained.  For  a 
detailed  mathematical  treatment  of  this  topic,  the  reader  is  referred  to 
Lee  [9]. 


CHAPTER  4 
APPLICATION  TO  AN  ADVERTISEMENT  PROBLEM 

In  this  chapter,  the  computational  aspects  of  this  technique  will  be 
discussed  with  respect  to  its  application  to  an  inventory  and  advertisement 
model  having  two  state  variables  and  one  control  variable. 

- 
4.1   DEVELOPMENT  OF  THE  MODEL 

The  diffusion  model  for  advertisement  was  originally  developed  by 
Teichroew  [14].   Consider  a  group  of  people  in  which  only  certain  members 
possess  a  particular  piece  of  information,  say,  about  a  manufacturing 
company's  product.   Suppose  that  the  total  number  of  persons  in  this  group 
remain  constant  and  that  the  diffusion  of  information  occurs  only  through 
personal  contact.  The  number  of  contacts  made  by  an  average  informed  person 
in  an  arbitrary  unit  of  time  is  given  by  a  contact  coefficient.   This  co- 
efficient is  same  for  all  members  of  the  group.   In  a  contact,  the  contactee 
receives  information  if  he  does  not  already  have  it;  if  he  already  has  it, 
the  contact  is  wasted  in  the  sense  that  it  did  not  increase  the  number  of 
informed  people. 

Let  Q(0)  =  Q  ■  number  of  informed  people  at  time  t. . 

N  =  total  number  of  persons 

c  =  contact  coefficient:  the  number  of  contacts  made  by 
c 

one  informed  person  per  unit  time. 
Q(t)  =  number  of  informed  persons  at  time  t. 
Q(t)/N  =  proportion  of  informed  persons  at  time  t. 
1  -  Q(t)/N  •  proportion  of  uninformed  persons  at  time  t. 
c  Q(t)  dt  =  contacts  made  during  a  time  interval  dt. 


The  increase  in  the  total  number  of  informed  people  during  a  short  interval 
of  time  At  is  obtained  by  multiplying  the  number  of  contacts  by  the  propor- 
tion of  uninformed  persons ,  because  an  increase  in  informed  members  is 
caused  only  by  contacts  with  uninformed  group.   Hence, 

dQCt)  -  c    Q(t)   dt    (1  -  Q(t)/N) 


49M«cc       QCt)   Cl-SM)  (25) 


Suppose  now  that  the  manufacturing  company  can  influence  the  number  of 
contacts  by  spending  money  for  advertising.   Specifically,  it  can  increase 
the  number  of  contacts  made  by  the  informed  people  by  an  additional  number 
A  per  unit  of  time.   Thus, 


4SM-QU)     Ccc  +  aco)     n-5!^-]  (26) 


If  each  informed  person  buys  c  units  of  the  company's  product  and 
if  SCt)  represents  the  sale  at  time  t,  then 


SCt)  -  cq   QCt)  (27) 


Let  c  =  1,  and  substitute  for  Q(t)  in  Equation; (26) , 


^f^1  =  S(t)    (cc  +  A(t))    [1  -  $&k  (28) 
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The  rate  of  change  of  the  company's  inventory,  I  CO,  is  given  by 


*$*■   =  P(t)  -  S(0  (29) 


where  P(t)  =  production  rate  at  time  t.  The  production  rate  is  assumed 
to  be  a  linear  function  given  by 

P(t)  =  a  +  bt 

where  a,  b  are  constants  and  t  is  time. 

This  is  a  typical  industrial  management  problem  where  the  management 
wishes  to  maximize  the  profit  given  by  Equation  (30) . 


J-  f  [c  SCO  -  cT  (I  -  I(t))2  -  cA  S(t)A2(t)]dt  (30) 

J  I   m  A 

i 
where  J  is  the  net  total  profit,  c  is  the  revenue  from  sale  of  one  unit  of 


The  role  of  the  management  in  this  particular  case  is  to  select  the 
optimal  policy  from  among  all  feasible  solutions  which  gives  the  maximum 
profit. 

4.2   DEFINITION  OF  THE  PROBLEM 

Maximize 

tf 
J  -|      [c  SCt)    -  Cl(Im  -  ICO)2  -  cA  S(t)    A2(t)]dt  (30) 
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subject   to 


F(t)   =  a  +  bt  (31) 

^r^-  P(t)    -   S(t)  (32) 

at 


^|~£i=   S(t)    (cc  +  A(t))    [1  -^~-]  (33) 


with  boundary   conditions 


I(t.)   -  1°  and  S(t   )   =   S°  (34) 


4.3   FORMULATION  OF  THE  PROBLEM 

The  above  optimization  problem  can  be  solved  by  calculus  of  variations 
with  the  help  of  the  quasilinearization  technique.   For  detailed  treatment 
of  the  calculus  of  variations,  the  reader  is  referred  to  Bliss  [3]  and 
Elsgolc  [5]. 

Equations  (31)  through  (34)  can  be  rewritten  as 


£^£  -  (a  +  bt)  +  s(t)  -  0  (35) 


^ir1  - s(t)  (cc  +  a^i1  -  ip-i  -  °  <36> 


We  have  two  state  variables,  I  and  S,  and  one  control  variable  A.   Intro- 
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duce  Lagrange  multipliers  X.,  X  and  constant  multipliers  8  ,  6„  and  define 
the  following  functions. 

Scc   SA 
F  =  [X    {1   -  (a  +  bt)  +  S)  +  A2(S  -  S(cc  +  A  -  — -)) 


+  cS  -  cT(I  -  I)2  -  c.  S  A2]  (38) 

1  m         A 


and 

G  »  [e1(i(0)  -  1°)  +  82(S(0)  -  S°)]  (39) 


where  the  notation  I  represents  the  first  differential  -j—  .   The  Euler- 
Lagrange  equations  [11] 


.1  <4Z)  _  _1Z  .  0  (40) 

dt  v3y.;   3y . 


and 


3A 


can  now  be  applied  to  Equation  (38)  to  obtain  relations  for  the.  Lagrange 


-^2(I(Im-D  (42) 


dA  2c  SX    2ASX 

^|=X1+c-cAA  -ccX2-X2A  +  -^  +  -F^  (43) 
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We  need  boundary  conditions  for  the  Lagrange  multipliers  X  ,  X«.  They 
were  obtained  by  applying  the  transversality  condition  [11]. 


3_G_ 


3E_ 


=  0     or 


ac 


it 

3y, 


=  o 


Applying  this  condition  to  Equations  (38)  and  (39) ,  we  obtain 


(44) 


Now  we  have  four  differential  equations  with  two  initial  and  two  final 
conditions,  which  make  the  problem,  a  two  point  boundary  value  type.   Applying 
condition  (41)  to  Equation  (38) ,  we  obtain 


A 


(43) 


Since  it  is  possible  to  express  the  control  variable  explicity  in  terms 
of  the  state  variables,  let  us  eliminate  A  from  all  the  performance  equations. 


dt 


a  +  bt  -  S 


2    2  3 

„         c  S*   S%  SX,   S% 

M  .   s  _  s   +    2 I ?- 

dt    c     N     c.N  2c.    ,  ..2 

A  A   2c .N 


(46) 


(47) 


-£•   2cT(I   -  I) 
at     1  m 


(48) 
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dX2       .   ,   +3S2X2+  X2     .   +2CCSX2 


(49) 


The  boundary  conditions  are  given  by  Equations  (34)  and  (44). 

4.4   QUASILINEARIZATION 

Observe  that  Equations  (47)  and  (49)  are  non-linear.  They  should  be 
linearized.   The  linearization  procedure  is  the  same  as  described  in 
Chapter  3.   Referring  to  Equation  (21),  we  need  the  expressions  for  f  ,  f,  , 


3g2    3g2 
This  matrix  can  be  obtained  from  Equations  (47)  and  (49) . 


2c  S   2SX„ 


*2   3S2X2 


2c, N     A 
A 


c.N   2c.    ,  „2 
A      A   2c,  N 

A 

4    ,  2ccA2 


c.N   2c.    ,  v2 
A      A   2c, N 
A 

3S2X.    \.        2SX.   2c  S 

-  c  +  T  +  o »t—  +  ~S — 

C   2c/  2cA   CAN     N 
A 


The  linearized  equations  and  the  recurrance  relations  can  be  developed  in 
accordance  with  Equations  (21)  through  (24) . 
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*ga  - .  +  bt  -  sn+1  (so) 


as  c  s  s  A„  s^\  s  X. 

n+1  _    ,  c  n         n  2,n         n  2,n  n  2,n   , 

~     lCoE„    _H+  r.    H  9^.  97~ia2"    J 


dt  l    c  n  N  cAN  2cA  2cANz 

2 
2c   s  2s  A_  A„  3s   A„ 

,    /  Nr  c  n    ,        n  2,n  2,n  n  2,n   , 

+   (sn+l  "  Sn}  [cc  "  ~T~  +  ~TJ-  '   2ct~  "   ?TW       ] 
AAA 

s2  s  s3 

+    ^2,n+1^2;n>tc7N-2-c7-rcV]  <51> 

A  A  A 


-^  -  2Vm  "  2cxVi  (52) 


dA  ,  3s   A  A  s   X,  2c  s  A. 

2,n+l        r,  ,  .        n  2,n    ,      2,n  n  2,n    ,        c  n  2,n   , 

j- —  -    [A,  +  c  -  c  A„       +  —j — =3—  +  -: — " *-  4 — *—  ] 

dt  l,n+l  c  2,n  4c, Nz  4c,  c.N  N 

AAA 

2  2 

3s  A.  A.  2c   A. 

,    /  u      D  2 ,n  2,ri  c  2,n   , 

+   (sn+l  "   8n)["2^i-  "  ^H~  +  ~~ T"^  ] 

2 

A  3s   A  2s   A.  2c   s 

+    (>2,n+l  "  X2,n)[2ct"  "   Cc  +  2cTn^ O-  +  T"  ]      (53) 

A  A  A 


The  boundary  conditions   are  given  by  equations    (34)    and  (44). 

Equations    (50)    through    (53)   are  ordinary   linear  differential  equa- 
tions and  with  the  boundary  conditions  given  by  equations    (34)   and    (44) , 
they  form  a  two  point  boundary  value  problem.      This  problem  can  now  be 
solved  by  superposition  principle.      Selecting   the  initial  conditions, 
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for  the  particular  and  homogeneous  solutions  such  that  they  satisfy  the 
given  initial  conditions,  the  general  solution  can  be  written  as 


I(t)  =  I  (t)  +  Ajl^t)  +  A2I2h(t) 

S(t)  =  S  (t)  +  AjSjjjCt)  +  A2S2h(t) 

Xj(t)  =  ^x^Ct)  +  *!*!  jjjW  +  A2xi  2h(t) 

X2(t)  =  A2)p(t)  +  AlA2,lh(t)  +  A2A2,2h(t)  (54) 

These  equations  are  derived  in  accordance  with  Equations  (11)  and  (12) . 

After  obtaining  the  final  solution  with  the  superposition  principle, 
Equations  (30)  and  (45)  can  be  solved  for  the  profit  and  advertisement, 
respectively.   This  completes  one  iteration.   Further  iteration  was  allowed 
until  desired  accuracy  was  obtained. 

4.5  NUMERICAL  ASPECTS 

In  order  to  solve  this  problem,  the  constants  were  assumed  to  have 

the  following  values . 

a  =  70  cA  =  1.5  1(0)  =  1°  =  20 

b  =  100  Cj  =  0.15  S(0)  =  S°  =  20 

cc  =  2  N  -  150  t.  -  0  t  =  1.0 

c  =  10  1-50  At  -  0.01 

m 

As  discussed  before,  we  need  initial  approximations  to  start  the 
solution.   Since  only  two  Equations,  (51)  and  (53),  were  non-linear,  we 
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obtained  from  intuition  and  knowledge  about  the  system.  Various  sets  of 
initial  approximations  used  in  this  problem  are  listed  in  Table  1. 

Solution  of  this  problem  by  the  superposition  principle  requires  a 
set  of  particular  solutions  and  four  sets  of  homogeneous  solutions.   How- 
ever, as  discussed  previously,  if  the  initial  values  for  the  particular 
and  homogeneous  solutions  are  chosen  such  that  they  satisfy  the  given 
initial  conditions,  only  two  sets  of  homogeneous  solutions  are  needed. 
The  following  set  of  values  at  the  initial  time  satisfy  this  condition 
and  hence  they  were  used  as  the  initial  values  for  the  particular  solution. 


I  (0)  =  20  X,   (0)  =  0.0 
P  1>P 

S  (0)  =  20  X,   (0)  =  0.0  (56) 

P  2,p 


The  initial  values  for  the  two  sets  of  homogeneous  solutions  were 
ssumed  as 


Iih(0> 

V0) 

Set  1 

0 

0 

Set  2 

0 

0 

xl,ih<°> 


1  0 

0  1    (57) 


4.6   COMPUTATIONAL  ASPECTS 

Using  the  initial  values  given  by  Equation  (56) ,  the  set  of  linear 
differential  Equations  (50)  through  (53)  were  solved  using  the  Runge-Kutta 
method  to  obtain  the  particular  solution. 

For  homogeneous  solutions,  the  known  terms  in  Equations  (50)  through 
(53)  were  set  to  zero  and  the  modified  equations  were  solved  by  Runge- 
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Table  1.  List  of  initial  approximations. 


Set  No.  SQ(t)  A2  Q(t) 


1 
2 
3 
4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 


450.0 

20.0 

350.0 

15.0 

300.0 

12.5 

300.0 

10.0 

275.0 

10.0 

250.0 

10.0 

200.0 

7.0 

150.0 

2.0 

100.0 

1.5 

50.0 

-  1.0 

28.0 

-  0.25 

25.0 

-  0.50 

22.0 

0.125 

20.0 

0.0 

5.0 

-  3.0 
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[READ  a,   b,    c,    c    ,    c, ,   H.C,,    I 
. _e  .    A^      '   I       m 


T 


Read  S.(e),  X2  Q(t) 


Eta 


PARTICULAR  SOLUTTON 

k^  =  2°  W0)  "  20 

Xl,n+1(0)  "  °  X2,n+1(0)  *  ° 

P  -  1       I  *  -i 


£ 


SOLVE  EQNS.    (51-54)    BiT 
SUBROUTINE  RKT 


W(t)  ■  WtJ 

Xlp,n+l(t)   =  W 
Vn+l(t)    -  X2,n+l(t) 


S2h,n+l(t)   *  W° 
X12h,n+l(t)   "   *l.«rU(t) 
X22K,n+l(t)   "    X2,n+l(t) 


HOMOGENEOUS  SOLUTION 


Xl,n+l(°> 

P  =  0 


1  V^M 


INTEGRATION  CONSTANTS 


55  using  Cramer's  Rule 


■*— 


I  >  0 


Jlh,n- 

Xllh,n+l(t>  "  Xl,,+1^> 

X21h,n+l(t)  *  X2,n+l(C) 


HOMOGENEOUS  SOLUTION 


"1 


xi.=u<°> 

?  -  0 


W0) " .° 

I_s_o   I 


FINAL  SOLUTION  USING 
EQUATIONS  55 


3— |n  =  n  +  l{~ 


Fig.  2.   Computer  Logic  Diagram  for  an  Advertisement  Problem. 
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Kutta  method  using  the  initial  conditions  given  by  Equation  (57).   These 
were  the  first  and  the  second  set  of  homogeneous  solutions.  Last  two 
Equations  in  (54)  at  final  time  t  ■  1  were  used  to  solve  for  the  two 
integration  constants,  A.  and  A  .   The  solution  was  obtained  by  Cramer's 
rule.  Next,  the  general  solutions  for  the  two  state  variables  and  two 
Lagrange  multipliers  were  obtained  by  using  the  superposition  principle, 
Equation  (54) . 

The  control  variable  A  and  objective  function,  J,  were  obtained  next 
using  Equations  (45)  and  (30)  respectively.  For  simplicity  the  following 
approximation  was  used  to  calculate  the  total  profit. 


f  2  , 

J  =  I      tcS(t)  -  cT(I  -  I(O)   -  cS(t)A/(t)]At 
t.  I  m  A 


In  general,  9  iterations  were  allowed. 

The  IBM  360/50  computer  system  was  used  for  all  these  computations. 
Computer  logic  diagram  is  shown  in  Fig.  2.  The  computer  program  is  given 
in  Appendix  2. 

4 . 7   RESULTS 

The  optimal  profit  in  this  program  was  J  ■  587.80  and  the  optimal 
initial  and  final  values  are 

1(0)  =  20  s(0)  =  20  A(0)  =  3.98 

1(1)  =  66.15         s(l)  =  115.69         A(l)  =  0 
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Out  of  the  15  different  sets  of  initial  approximations  listed  in  Table  1, 
the  first  five  sets  did  not  produce  convergence. 

In  set  1,  the  particular  and  homogeneous  solutions  of  all  the  four 

40 

variables  at  the  final  time,  t  ,  involve  terms  of  the  order  10   .  As  a 

result,  the  calculation  of  integration  constants  by  Cramer's  rule  involves 

80 
terms  of  the  order  10   ,  which  cannot  be  handled  by  IBM  360  computer  and 

exponential  overflow  was  resulted. 

Set  2  encountered  a  similar  problem.   The  lowest  term  in  the  partic- 

13 
ular  and  homogeneous  solutions  at  final  time  t,,  was  of  the  order  10 

This  did  not  cause  any  difficulty  in  the  calculation  of  integration  con- 
stants, but  resulted  in  exponential  overflow  in  the  computation  of  final 
solution  of  first  iteration. 

Sets  3,  4,  and  5  encountered  basically  the  same  problem.   For 
explanation,  results  of  set  4  are  used  here.  In  the  final  solution  of 
first  iteration,  the  following  results  were  obtained. 


A(l)  =  -0.187  x  1014      1(1)  =  0.706  x  107 
s(l)  =  -0.141  x  109 


as  a  result  of  such  large  numbers ,  the  computer  experienced  exponential 
overflow  and  stopped  computation  while  calculating  the  particular  solution 
of  the  second  iteration. 

Sets  6  through  15  converged  to  the  same  optimal  solution  in  about 
4-5  iterations.   The  convergence  rates  of  sales,  inventory,  and  advertisement 
are  shown  in  Figs.  3  through  11.  The  initial  approximations  used  are  sets 
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10,  12,  and  14  in  Table  1.  Fig.  12  shows  the  convergence  rate  of  the 
profit  function  with  set  14  of  Table  1  as  the  initial  approximation.  The 
convergence  rates  of  the  initial  and  final  values  of  the  variables  are 
tabularized  in  Tables  2  through  5.   The  IBM  360/50  computer  took  about 
3.72  minutes  to  complete  9  iterations  for  this  problem  with  WATFOR  compiler. 

4.8   DISCUSSION 

The  results  show  that  this  problem  converged  with  ten  different 
and  far  from  optimal  initial  approximations.  The  optimal  curves  show  that 
the  profiles  of  the  state  and  control  variables  were  either  monotonically 
decreasing  or  monotonically  increasing.   This  made  qussilinearization 
method  more  effective. 

It  was  observed  from  Tables  2-5  that  convergence  was  obtained  in  4-5 
iterations  for  all  sets  of  initial  approximations  that  converged.  It  was 
concluded  that 

1.   The  quasilinearization  method  converges  quadratically, 
whenever  it  converges, 
and  2.  The  convergence  rate  is  almost  independent  of  the  choice 
of  initial  approximation,  if  the  latter  values  are 
within  the  convergence  interval  or  range. 
A  note  on  the  choice  of  initial  approximation  is  in  order.   In  this 
problem,  the  optimal  solution  of  sales  is  between  20.0  and  115.69.   But 
any  initial  approximation  of  sales  between  5.0  and  250.0  would  converge 
to  the  optimal  solution.   Hence,  choosing  this  value  should  not  be  a 
problem.   The  author  feels  that,  in  general,  the  basic  knowledge  of  physical 
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Fig.  3.   Convergence  Rate  of  Sales  in  an  Advertisement  Problem 
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Fig.  6.   Convergence  Rate  of  Inventory  in  an  Advertisement  Problem. 
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Fig.  7.   Convergence  Rate  of  Inventory  in  an  Advertisement  Problem 
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Fig.  8.   Convergence  Rate  of  Inventory  in  an  Advertisement  Problem. 
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Fig.  9.   Convergence  Rate  of  Advertisement  in  an  Advertisement  Problem. 
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Fig.  10.   Convergence  Rate  of  Advertisement  in  an  Advertisement 
Problem. 
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Fig.  11.   Convergence  Rate  of  Advertisement  in  an  Advertisement 
Problem. 
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behavior  of  the  system  is  enough  to  make  correct  choice.   Furthermore,  many 
numerical  schemes  have  been  devised  to  overcome  the  convergence  problem. 
One  such  scheme  is  the  data  perturbation  technique  [4]. 

In  order  to  further  investigate  the  convergence  and  other  computational 
aspects  of  this  problem,  the  following  constants  were  used. 


a .=  0.7  cA  =  1.0  1(0)  =  1°  =  0.2 

b  -  1.0  c,  ■  0.15  S(0)  -  S°  =  0.2 

c  -  2.0  S  -  1.5 

c 

c  =  10.0  I  ■  1.0  it  =  0.01      (58) 

m 


The  initial  approximations  used  were 


The  initial  values  for  the  one  particular  and  two  homogeneous  solutions 
were 


K0) 

8(0) 

Xx(0) 

x2(o: 

Particular  soln. 

0.2 

0.2 

0.0 

0.0 

Homo.  soln.  set  1 

0.0 

0.0 

1.0 

0.0 

Homo.  soln.  set  2 

0.0 

0.0 

0.0 

1.0 

The  convergence  rates  are  shown  in  Table  6.   The  problem  converged  In  4 
iterations.   The  optimal  profiles  of  I,  S,  and  A  are  shown  in  Fig.  13. 
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Table  6.   Convergence  rates  for  the  modified  problem,  Equation  (58) 

Iteration  Time  }     (t)     (t)     A(t)     j 

Number     t  •  1     '  2 


0.0 

0.20 

0.0 

0.25 

0.20 

0.0 

0.50 

0.20 

0.0 

0.75 

0.20 

0.0 

1.00 

0.20 

0.0 

0.0 

0.200 

0.200 

-0.225 

-0.225 

9.763 

0.25 

0.292 

0.714 

-0.169 

-13.563 

3.553 

0.50 

0.313 

1.284 

-0.118 

-7.337 

0.529 

0.75 

0.240 

1.985 

-0.064 

-3.010 

-0.487 

1.00 

0.027 

2.915 

-0.0 

0.0 

0.0 

0.0 

0.200 

0.200 

-0.187 

-13.903 

6.025 

0.25 

0.308 

0.592 

-0.131 

-6.851 

2.072 

0.50 

0.383 

0.945 

-0.082 

-3.526 

0.653 

0.75 

0.448 

1.146 

-0.038 

-1.985 

0.234 

1.00 

0.532 

1.433 

0.0 

0.0 

0.0 

0.0 

0.200 

0.200 

-0.180 

-12.369 

5.360 

0.25 

0.313 

0.552 

-0.124 

-6.308 

1.994 

0.50 

0.403 

0.861 

-0.076 

-3.994 

0.851 

0.75 

0.490 

1.078 

-0.034 

-2.173 

0.306 

1.00 

0.595 

1.223 

0.0 

0.0 

0.0 

0.0 

0.200 

0.200 

-0.180 

-12.411 

5.378 

0.25 

0.313 

0.552 

-0.124 

-6.333 

2.001 

0.50 

0.403 

0.860 

-0.076 

-4.015 

0.856 

0.75 

0.490 

1.077 

-0.034 

-2.175 

0.307 

1.00 

0.596 

1.219 

0.0 

0.0 

0.0 

0.0 

0.200 

0.200 

-0.180 

-12.410 

5.378 

0.25 

0.313 

0.552 

-0.124 

-6.332 

2.001 

0.50 

0.403 

0.860 

-0.076 

-4.015 

0.856 

0.75 

0.490 

1.077 

-0.034 

-2.175 

0.307 

1.00 

0.596 

1.219 

0.0 

0.0 

0.0 

8.631 


6.682 


6.669 


6.669 
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Fig.  13.   Optimal  Profiles  of  I(t),  S(t)  and  A(t)  for  modified 
Problem,  Eqn.  (59). 
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Most  important  advantage  of  this  technique  is  that  the  control 
variable  can  be  eliminated  from  the  performance  equations.   Hence,  initial 
guess  for  the  control  variable  is  not  required  in  solving  for  the  state 
variables.   For  this  reason,  this  technique  is  superior  to  the  gradient 
technique  because  any  error  in  guessing  the  initial  control  can  result 
in  failure  to  obtain  the  solution. 

It  may  be  interesting  to  note  that  the  results  of  this  problem, 
using  the  initial  values  listed  in  Equation  (55),  compare  favorably  with 
the  results  obtained  by  the  first  variational  gradient  technique  [15]. 
The  modified  problem  with  numerical  values  given  by  Equation  (58)  has 
also  been  solved  by  the  second  variation  technique  [13] .   The  present 
results  again  compare  favorably  with  this  result. 

It  will  be  observed  in  the  next  chapter,  that  even  if  the  control 
variable  cannot  be  canceled  from  the  performance  equations,  quasilineari- 
zation  method  still  works  well. 
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CHAPTER  5 
APPLICATION  TO  AN  ADVERTISEMENT  AND  PRODUCTION  PROBLEM 

We  now  wish  to  apply  the  quasilinearization  technique  to  a  more  com- 
plex problem,  namely  an  advertisement  and  production  problem.   This  problem 
has  six  state  variables  and  three  control  variables.  In  addition,  the 
profiles  are  fairly  unstable  due  to  the  rapid  change  of  the  variables 
with  time. 

5.1   DEVELOPMENT  OF  THE  MODEL 

Consider  the  manufacturing  process  shown  in  Fig.  14.  There  are  two 
chemical  reactors  in  which  the  following  consecutive  reactions  take  place 

A >  B  >  C 

Both  these  reactions  are  first  order.   The  component  B  is  the  desired 
product  and  C  is  the  waste  product.   Suppose  B  is  a  new  product  which  needs 
advertisement  to  boost  the  sales.   Furthermore,  to  protect  against  fluct- 
uations in  demand,  an  inventory  will  be  assumed  for  B.  A  and  C  are 
assumed  to  have  unlimited  market  at  fixed  price  and  they  are  sold  as  soon 
as  manufactured. 
Let  x  and 
respectively.   Under  steady  state  conditions,  from  material  balance,  we  have 


amount  of  _  amount  of   amount  of  A 

A  in       A  out     transformed  to  B  ^     ' 

amount  of  _  amount  of   amount  of  B      _  amount  of  B 

B  in        B  out     transformed  to  C  "  produced  from  A     ^60^ 
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v.  "  volume  of  chemical  reactor  i,  i  =  1,  2. 

q  =  flow  rate 

k  .  =  reaction  rate  constant  of  the  first  reaction  in  reactor  i 
ai 

Vl  .   =   reaction  rate  constant  of  the  second  reaction  in  reactor  i 

G  ,  G,  =  frequency  constants  of  the  first  and  second  reactions,  respectively 
a   b 

E  ,  E,  ■  activation  energies  of  the  first  and  second  reactions,  respectively 
a   d 

R  =  gas  constant 


The  kinetics  of  the  reactions  can  now  be  written  as 


q(X(3  -  Xj)    -   V^^  =  0 


at  steady  state.   Under  unsteady  state  conditions  we  have 


dx. 

vl  dT  =  q(x0  '  xl)  "  VlkalXl  C61) 


similarly,  eqn.  (60)  can  be  rewritten  as 


q(y0  -  yx)  -  v^^  +  v1kalx1  =  o 
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at  steady  state.  Under  steady  state  conditions,  we  have 

dyl 
vl  IT  =  q(y0  -  yl}  "  vlkblyl  +  vlkalxl  (62) 

With  similar  arguments,  the  kinetics  of  the  reactions  in  the  second 
reactor  can  be  written  as 


dx2 
V2  dt~=  q(xl  "  X2}  _  V2ka2X2  (63) 

and 


V2  IT  =  q(yl  -  y2)   -   V2kb2y2  +  V2ka2X2  (6A) 


The  reaction  rate  constants  are  defined  as 


al    a 


Ea  *b 


\i m  Ga ex? (- 1& '     \2'%  exp <-  SrJ  <«> 


As  indicated  before  B  is  the  desired  product  and  it  needs  inventory 
and  advertisement.   The  performance  equation  for  the  inventory  is 


rate  of  change  _  production   sales 
of  inventory       rate      rate 


£•<&■>-   ■  (66) 
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where  I  is  the  inventory  and  S  is  sales. 

The  performance  equation  for  advertisement  is  the  same  as  Equation 
(26)  in  the  last  chapter.   Again,  let  c  »  1.  According  to  Equation  (27), 
we  have 


Thus, 


S(t)  ■=  c  Q(t)  =  Q(t) 


af=S(cc  +  A)[l-f]  <67> 


where  c  is  the  contact  coefficient.  A  is  the  advertisement  and  N  repre- 
sents the  total  number  of  people  in  the  group. 

Equations  (61),  (62),  (63),  (64),  (66),  and  (67)  describe  the  system 
letely.   We  have  six  state  variables 
three  control  variables,  T  ,  T„,  and  A. 

The  management  in  this  particular  industrial  system  is  confronted  with 
the  problem  of  selecting  three  control  variables  such  that  the  following 
profit  function,  J,  is  maximized. 

profit  =    [revenue  of  B  +  revenue  of  A  +  revenue  of  C 
Jti 

inventory   advertisement   manufacturing, 
cost  cost  cost 

Mathematically, 


-Jt  [clCqS+C2" 


,q*2  +  c3q(l  -  x2  -  y2)  -  Cl(Im  -  I)2 


i 

.2  2 


5.2   DEFINITION  OF  THE  PROBLEM 
Maximize  the  functional 


J  "  J   IclcqS  +  c2^2  +  c3q(1  -  x2  "  V  -  cI(Im  "  T) 
i 
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CAA"S"  "  CT{(Tln.  "  V  +  (T1  "  V  }]dt  (68) 


"  ^^  "  CT{(Tlm  "  Tl)2  +  (X1  "  T2)2}]dt  (69) 

subject  to  the  constraints  of 

a* 

VldT=  q(x0  -  V   -  vikalxi  W) 

dy2 
Vl  W  "  *%  -  7j)   -  v1ktly1  +  VlkalXl  (71) 

dx2  . 

V2  di~  =  qK7-l  '  V    *  V2ka2X2  <72> 

dy2 

V2  dT  =  q(yI  "  y2}  '  v2kiay2  +  V2ka2x2  <73> 

dl 

dt  =  qy2  "  S  (?*) 


||=S(cc  +  A)[l-|]  (75) 


with  boundary  conditions 


59 


X1(V  =  Xl         W  =  y2 

y/t.)  =  y°         I(t.)  -  1°  I(tf)  =  I1 

x2(t.)  -  x°         Stt±>  =  S°  (76) 


5.3   FORMULATION  OP  THE  PROBLEM 

It  was  required  to  find  the  optimal  value  of  the  state  variables  and 
control  variables  so  that  the  objective  function  is  maximized.   This 
problem  can  be  solved  by  calculus  of  variations.  The  procedure  for 
obtaining  the  solution  remains  essentially  the  same. 

Equations  (70)  through  (76)  can  be  rewritten  as 


J.  -  a~  (xA  -  x,)  +  G  e  RT,"  x  =  0  (77) 

1   v.   0    1     a      1   1 


h  -   vx  (y0  -  yl>  +  Gb  e  RT1  yl  ~  Ga  e  RT1  Xl  "  °  (78) 


X2  -  v^  (X1  -  X2)  +  Ga  e"  RT2  X2  =  °  (7S) 

Eb  Ea 

y2  "  v^    (yl  "  y2>  +%e'W2     y2  "  Ga  e"  ^     x2  "  °  (80) 

(81) 


S  -  (ccS  +  AS)[1  -  |]  =  0  (82) 

dx. 


The  symbol  x  represents  -r— 
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Introduce  lagrange  multipliers,  ^ . ,  i  "  1,  .  .  . ,  6,  and  constant 
multipliers  9  ,  j  =  1,  .  .  .,  7,  and  define  the  following  functions. 

E 
F  -    [Al(il  -  4-   (x0  -  xx)  +  Ga  e"  ly  XX) 

E  E 

*  V*i  -  *-  (y0  -  yj)  +  S  e"  ^  yi  "  Ga  e"  ^   V 


+  S(i2  -  ^   (X1  -  V  +  Ga  e"  RT2     X2 

Eb  E 

+  \(y2  -  ^  <y1  -  y2>  +  %  ."  ht^  y2  -  ca  r  ^  ^ 

+  X,(j  -  qy,  +  S) 


+  A6(S  _  cS  -  AS  +  ^f-  +  A|_> 


+  CLCqS+  C2qX2  +  C3q(1  "  X2  "  V   "  C!(Im  "   I)2 


-  cAA2S2  -  ^{(T^  -  T^2  +    IJ1  -  T2)2}]  (83) 


g  -  [ejfejfcj)  -  x°)  +  e2(yi(t±)  -  y°)  +  e3(X2(t.)  -  x°) 


+  VW  -  y°)  +  95<i(tt)  -  i°) 


+  e6(i(tf)  -  i1)  +  e7(s(t.)  -  s0)]  (84) 


The  Euler  -  Lagrange  equations  [  11  ] , 
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d   ,3F  .    3F    _  ,„% 


and 

3F 


3Z  "°  -x  <86> 


can  now  be  applied  to  Equations  (83)  to  obtain  the  relationships  for  the 
six  lagrange  multiplier  equations. 


dX      A    A  _  _b 

F-^tJ*^'   rti  .     (88) 


<U-   X„q  a 


dX    X  q  % 

IT  m  T~  +  Y'b  e"  RT2  -«««3  +  V  (90) 


dX 

^-  =  2cTI  -  2cTI  (91) 

dt      I  m     I 


dXfi  2cSX    2ASX,       - 

3T  ■  Cl  +  A5  "  CA6  -  W6  +  -IT1  +  -T1  -  2CAA  S  <92> 


Application  of  Equations  (86)  to  (83)  yields 
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A  =  —  f  —  -  —  )  (93) 

A   2c,   N   S  *  **  ' 

A 

(X.  -  AjG  E  x     \  l.GE,     \_ 

_1 2  a  a  1  e  +  ^  b  °  X  e  RT   -  2c  (2T  -T  -T  )  =  0    (94) 

RT  RT1 

! 

(X„  -  X.)G  E  x,    \_  WKy,  _  ^b_ 

—3 Y^-5-^  e~  RT2  +  *  b  °  2  e  RT2  +  2cT(T1-T2)  =  0       (95) 


Equation  (93)  gives  explicit  expression  of  the  control  variable  A. 
Hence  A  can  be  eliminated  in  all  the  performance  equations.   However,  the 
control  variables,  T1 
cannot  be  eliminated. 

Substituting  the  expressions  for  A  into  Equations  (82)  and  (92),  we 
obtain 

.2   SA,   X.  S2\, 

||  -  CS  -  Sf  +  -|  -  T6- 6-j  (96) 

dt        N    CAN   2CA   2c  N2 
A 

dX.  ll         2cSX,   sx, 

A  2c.  N 

r'  A 

Notice  that  these  two  equations  are  non-linear. 

Equations  (77)  through  (81),  (87)  through  (91),  (96)  and  (97)  represent 
the  system.   For  these  12  differential  equations  we  have  only  7  boundary 
conditions  given  by  Equation  (76).  The  additional  5  boundary  conditions 
can  be  obtained  by  applying  the  transversality  condition  [11] . 
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3G 

3y, 


tf< 


=   0 


i    It, 


3G 

¥y. 


3F 


=   0 


to  Equations    (83)    and    (84) 


\(t£)  -  0 

W  -  0 

»,<V  -  0 

W  =  o 

w "  ° 


(98) 


The  boundary  conditions  given  by  Equations  (76)  and  (98)  make  this 
system,  a  two  point  boundary  value  problem. 

Let  us  consider  the  case  that  the  final  condition  on  the  inventory 
was  not  given.  Equations  (69)  through  (75)  remain  unchanged.  Equation 
(76)  can  be  modified  as 


w  ■  *i 


K2(t.) 


y2(tt)  ■  y2 


s(tt)  =  su 


(99) 


Definition  of  the  function  F,  Equation  (83),  remains  the  same,  but  the 
function  G  is  modified  as 
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g  =  [e1(x1(t.)  -  x°)  +  e2(y1(t.)  -  y°)  +  s^C*^  "  "a5 


+  e4(y2(t.)  -  y°)  +  e5(i(tl)  -  i°) 


(100) 


Since  E  remains  unchanged,  Equations  (87)  through  (97)  are  valid  in 
this  case. 

Applying  the  transversality  condition  [11]  to  Equation  (100),  we  have 


X1(tf)  =  0 


0  -  A4(tf)  -  0 

0  -  x5(tf)  =  0 


(103) 


5.4   QUASILINEARIZATION 

Only  Equations  (96)  and  (97)  are  non-linear.  The  linearization 

procedure  is  the  same  as  described  in  Chapter  3.  Referring  to  Equation 

(21),  we  need  the  expressions  for  f  ,  f,  ,  g  ,  g,  .  In  other  words,  we  need 


»g,    3g, 


3g2   382 
3S~   SU 
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This  matrix  can  b 

;  obtained  from  Equations  (96)  and  (97). 

2cS 

.  X6   SA6       SI      S2 
cAN   ^N?  •     cAN   2cA  "  2c^P 

J  = 

2cA, 

2 
X6               2cS      X6     X6 

_ 

(102) 

!c,N2         '      N    C   c.N   c.N^ 
A                         A     A 

Linearization  and 

recurrence  relations  were  developed  in 

accordance  with 

equations  (21)  through  (24) . 

l,n+l   q 

E 
a 

(x0  "  Xl,^   -   Ga  e"  RT1  Xl,n+1 

(103) 

dt     v 

dyl,n+l   q_ 
dt     v. 

(y0  "  yi.n+l5  "  Gb  6  RT1  yl,n+l  +  Ga 

E 
*~A    xl,n+l  W) 

2,n+l   q 
dt      v2 

E 
a 

(xl,n+l  "  *2,a«>  "  Ga  e_  RT2  X2,n+1 

(105) 

dy2,n+l   q 

Eb 
(yl,n+l  "  *2,n«>  "  Gb  e  RT2  y2,n+l  + 

E 
Ga  e"  R^  X2,n+1  (106> 

dt      v2 

dIn+l 

~1T  =  qy2,n+l  _  Sn+1 

(107) 

-dT  =  tcSn 

2                  2 
cS    S  X,     X,     SI, 
n    n  o,n    6,n    n  6,n  -, 

N      c.N    2c.   "  ,  „2   J 
A       A    2c, N 
A 

+  (Sn+l 

„  ,,  r    2cSn  ,   6,n    n'  6,n 

VlC    N     cAN  "  c„N2 
A       A 
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S  s2 

+    &fi  -.1    -  ^  J^  -  ^Tr S"J  <108> 

o ,  nrl  o ,  n     c.N        ^c.        „     vtz 


A  'A       2c.  N 

A 


E 


l,n+l  ,   l.n+1         3,n+l*    ,    ,,  ,  >   _       -  ==- 


5>n+1  =   2c   I     -   2c   I 
dt  I  ■  I  n+1 


dX,      .,                                                           \\  2cS  X,  S  A? 

6, n+1        r         .    ,                      ,                 6,n    .          n  6,n  n  6,n   , 

dt                1          5,n+l            6,n       2c  J              N  2 
2cA,            X, 


(109) 


^  ■  '^  -  ^f  >  ♦  A2>n+1Gb   ."  B"  (110) 

E 
_3jEtl  .  _ABil  +  (X3>n+i  .  X4)n+i)Ga  e-  RT2    +  q(c2  _  ^  (111) 


4, n+1  4, n+1   .   ,  _  :r=—  ,  /ins 

—dt-  =  — ^-  +  X4,n+lGb  e     EI2     "  c3q  "  «X5,*U  (112) 


n+1    n'1   N      ,  ..2 
2c. N 
A 


2cS        X,     S  X, 

+  ft.   -  -X,      )[-5-i"  «--V  +  -J^l*S]  (I") 

6, n+1    6,n   N        c.N    c.N-' 
A      A 


The  boundary  conditions  are  given  by  Equations  (76)  and  (98)  or  (76)  and 
(101). 

Equations  (103)  through  (113)  are  ordinary  linear  differential  equa- 
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tions  and  with  the  boundary  conditions  given  by  equations  (76)  and  (98) 
or  (76)  and  (101),  they  form  a  two  point  boundary  value  problem.   This 
problem  can  now  be  solved  by  the  superposition  principle.   If  the  initial 
values  for  the  particular  and  homogeneous  solutions  are  selected  such  that 
they  satisfy  the  initial  conditions,  the  general  solution  can  be  given 
as 


(114) 


Xl(t)  =  x   (t)  +  IVl,kh(t) 
r  k=l 

yi(t)  -  y  (t)  +  I  Vl,kh<°  U15) 

r      k=l 


(116) 


x2(t)  =  x2p(t)  +  k^V2,kh(t) 

y2(t)  -  ,  (t)  +  I  V2,kh(t)  (117) 

r      k=l 


I(t)  -  I  (t)  +  I   Vkh(t) 
r     k=l 

6 
S(t)  -  S  (t)  +  I  Vkh(t) 
r     k=l 

6 


(118) 


(119) 


X.(t)  -  Xip(t)  +  ^Vi.fcfcW  t  -  1.  ....  6.     (120)-(125) 

These  equations  are  derived  in  accordance  with  Equations  (11)  and  (12) . 

After  obtaining  the  solution  for  the  6  state  variables  and  6  Lagrange 
multipliers  by  the  superposition  principle,  Equations  (69),  (93),  (94),  and 
(95)  can  be  solved  for  the  profit  and  the  three  control  variables,  A,  T., 
T. ,  respectively. 
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All  these  calculations  complete  one  iteration.   Further  iterations 
were  allowed  until  desired  accuracy  was  achieved. 

5.5  NUMERICAL  ASPECTS 

Depending  upon  the  value  of  the  constants  and  the  boundary  conditions, 
this  problem  was  divided  into  classes  A,  B,  C,  D,  and  E.   The  object 
was  to  investigate  the  convergence  and  other  computational  aspects  of  this 
technique  from  different  angles. 

Problem  A 

The  following  values  were  assumed  for  the  various  parameters 


Ga  =  0.535  x  10   per  minute        N  «  100 

18 

G^  =  0.461  x  10   per  minute        c  =  1 

E  =  18000  cal/mole 
E.  =  30000  cal/mole 
R  ■  2  cal/mole  °K 
q  =  60  gal/min 


10  gallons  x.(t)  =  0.53    t  ■  0.0 

U  l 

it  -  0.01         yQ(t)  =  0.43    tf  =  1.0      (126) 

The  boundary  conditions  were 

x1(0)  =  0.53   yi(0)  -  0.43   x2(0)  =  0.53   y2(0)  -  0.43  1(0)  =  1.0 

S(0)  =  0.1     1(1)  =  10.0 

It  should  be  emphasized  that  class  A  was  the  only  problem  which  had 


69 


final  condition  on  the  inventory. 

There  were  only  two  non-linear  differential  equations,  hence  only 
two  initial  approximations,  SQ(t)  and  X,  n^  '   were  required.  The 
equations  for  the  control  variables,  T  and  T„,  are  implicit  and  cannot  be 
solved  directly.   Hence,  the  initial  approximations  for  T  and  T„  were 
also  required.   The  various  sets  of  initial  approximations  used  for  problem 
A  are  listed  in  Table  7. 

Problem  B 

The  same  parameters  used  in  problem  A  were  used  here,  except  that  the 
final  condition  on  the  inventory  was  removed.  As  a  result  of  this  change, 
according  to  Equation  (101),  X-d)  «  0.  The  boundary  conditions  are 

x  (0)  -  0.53  y,(0)  ■  0.43  x2(0)  =  0.53  y,,(0)  =  0.43  1(0)  >  1.0  S(0)  =  0.1 

As  in  the  last  case,  only  SQ(t),  Xfi  Q(t),  T±   Q(t),  and  T2  Q(t)  were  required 
as  the  initial  approximations.   They  are  listed  in  Table  8. 

Problem  C 

Some  of  the  parameters  were  changed.   For  clear  understanding  all  of 
them  are  rewritten  in  the  following 

G  ■  0.535  x  1011  per  minute        N  =  100 
a 

18 
G,  =  0.461  x  10   per  minute       c  =  1 

E  =  18000  cal/mole  (^  =  0.0005  $/°k 

BL  =  30000  cal/mole  cA  =  0.0002  $ 
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R  =  2  cal/mole  °k  c  =  5.0  $ 

q  =  60  gal./min.  c_  =  c-  ■ 

v  =  v.   -   12  gallons  c  =1.0 

I  =  20  gallons  x_(t)  =  0.53   t.  =  0.0 


The  boundary  conditions  were 

Xl(0)  =  0.53  yi(0)  =  0.43  x2(0)  =  0.53  y2(0)  =  0.43  1(0)  =  8.0  S(0)  =  0.1 

A  list  of  initial  approximations  is  shown  in  Table  9. 

Problem  D 

c,  =  0.01  S(0)  =  1.0 

A 

All  other  parameters  were  the  same  as  in  problem  C.   Three  different 
initial  approximations  were  used  for  this  problem.   They  are  tabulated  in 
Table  10. 

Problem  E 

The  only  difference  between  problems  E  and  C  is  in  the  initial  condition 
of  the  sales.   In  the  present  problem,  the  initial  condition  for  sales  is 

S(0)  =0.1 

All  others  values  are  the  same  as  in  problem  C.   The  values  of  the  initial 
approximation  used  are  given  in  Table  11. 

The  initial  values  used  for  the  particular  and  homogeneous  solutions  are 
given  in  Table  12. 
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Table  7. 

Initial  approximations 

for  problem  A. 

Set  No. 

Ti,o(t> 

T2,0(t) 

SQ(t) 

Vo(t> 

1A 

345 

345 

1 

0 

2A 

345 

345 

50 

-0.5 

3A 

345 

345 

30 

-1.0 

4A 

330 

330 

1 

0 

Table  8. 

Initial  approximations 

for  problem  B. 

Set  No. 

Tl,0(t> 

T2,0(t) 

yo 

X6,0(t) 

IB 

330 

330 

i 

0 

2B 

340 

340 

i 

0 

3B 

345 

345 

i 

0 

Table  9. 

Initial  approximations 

for  problem  C. 

Set  No. 

Tl,0(t> 

T2,0(t) 

SQ(t) 

X6,0(t> 

1C 

330 

330 

1 

0 

2C 

345 

345 

1 

0 

Table  10. 

Initial  approximations 

for  problem  D. 

Set  No. 

Tl,0(t> 

T2,0(t) 

s0(t) 

Vo(t> 

ID 

355 

355 

20 

-0.5 

2D 

345 

345 

50 

0 

3D 

350 

350 

25 

-0.25 

Table  11. 

Initial  approximations 

for  problem  E. 

Set  No. 

Ti,o(t) 

T2,0(t)  . 

yo 

X6,0(t) 

IE 

345 

345 

50 

0 
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Table  12.   Initial  conditions  used  for  obtaining  particular 
and  homogeneous  solutions. 


P.I. 

for 
A,  B 

P.I. 

for 
C,  E 

P.I. 

for 
D 

Homogeneous 

Solutions 

Set  1 

Set  2 

Set  3 

Set  4 

Set  5 

Set  6 

Xl(0) 

0.53 

0.53 

0.53 

0 

0 

0 

0 

0 

0 

7l(0) 

0.43 

0.43 

0.43 

0 

0 

0 

0 

0 

0 

x2(0) 

0.53 

0.53 

0.53 

0 

0 

0 

0 

0 

0 

y2(0) 

0.43 

0.43 

0.43 

0 

0 

0 

0 

0 

0 

KO) 

1.0 

8.0 

8.0 

0 

0 

0 

0 

0 

0 

S(0) 

0.1 

0.1 

1.0 

0 

0 

0 

0 

0 

0 

1<« 

0 

0 

0 

1 

0 

0 

0 

0 

0 

2(0) 

0 

0 

0 

0 

1 

0 

0 

0 

0 

3(0) 

0 

0 

0 

0 

0 

1 

0 

0 

0 

4(0) 

0 

0 

0 

0 

0 

0 

1 

0 

0 

5(o) 

0 

0 

0 

0 

0 

0 

0 

1 

0 

««» 

0 

0 

0 

0 

0 

0 

0 

0 

1 
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5.6   COMPUTATIONAL  ASPECTS 

Basically  the  same  procedure  was  followed  for  all  the  five  problems. 
This  procedure  was  essentially  the  same  as  that  used  in  Chapter  4.   With 
the  initial  values  given  in  Table  12,  a  set  of  particular  solutions  and 
six  sets  of  homogeneous  solutions  were  obtained  by  numerical  integration 
using  the  Runge-Kutta  method. 

In  order  to  solve  for  six  integration  constants,  six  equations  for 
which  the  final  conditions  were  known,  were  selected  from  Equations  (114) 
through  (125).   To  solve  this  6x6  matrix  on  computer,  matrix  inversion 
subroutine  SIMQ  supplied  by  IBM  was  used.  A  printout  of  this  subroutine 
is  shown  in  Appendix  3.  Using  these  integration  constants  with  the  newly 
obtained  particular  and  homogeneous  solutions,  the  final  solutions  for 
all  twelve  variables  were  obtained.   Next,  using  Equation  (93),  values  of 
advertisement  at  all  grid  points  were  obtained. 

For  simplicity  the  following  approximation  was  used  to  calculate  the 
total  profit. 

J  =  I    [CjC  S  +  c2qx2  +  c3q(l  "  x2  -  y2)  -  (^(1,  -  I) 


q 

"i 


cAA2S2  -  cT{(Tlm  -  Tx)2  +  dj  -  T2)2}]At  (128) 


Computation  of  T.  and  T„  were  rather  difficult  because  Equations  (94) 
and  (95)  are  implicit  in  T  and  T„.   To  overcome  this  difficulty,  c  was 
assumed  to  be  zero  and  the  last  term  in  both  the  equations  dropped  out. 
As  a  result,  explicit  expressions  were  derived  as 
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(E  -  E,  )/R 
T,  =   a    b 


1      G  x.  E  (X  -  A.) 
ln[-a  1  a  2 M 


(E  -  E,  )/R 
T  -    *    b 


2      G  x,  E   (X,  -  X,) 

^i  2  v  « \     ]  (l29) 

Gb  y2  Eb  X4 


These  equations  can  be  solved  easily,  but  examining  these  equations  care 
tf ,  the  denominator  would  involve  lnbr 


fully,  at  final  time  t  ■  tf,  the  denominator  would  involve  lnbd  which  is 


indeterminate. 

Another  approach  to  this  difficulty  was  to  apply  the  Newton-Raphson 
method  of  root  finding.   This  method  is  described  in  Appendix  1.  As 
indicated  earlier,  it  is  fundamentally  the  same  as  quasilinearization. 
To  start  the  Newton-Raphson  method,  initial  approximations  for  both  T  and 
T,  were  assumed  to  be  350.0.   With  these  values,  T  and  1     at  the  first 
grid  points  of  the  first  iteration  were  solved.  At  other  grid  points, 
solutions  of  T  and  T.  at  previous  grid  points  of  same  iteration  were  used 
as  the  initial  approximations.   Same  procedure  was  repeated  for  the  next 
iteration. 

This  scheme  encountered  convergence  problem  in  problems  3A,  3B,  ID 
and  IE.   The  logic  was  slightly  modified  to  overcome  this  trouble.  The 
procedure  remained  the  same  for  the  first  iteration.   For  other  iterations, 
solutions  of  T.  and  T  at  the  same  grid  point  of  previous  iteration  were 
used  as  the  initial  approximation.   With  this  change  the  convergence 
difficulty  was  overcome  in  problems  3B  and  IE.   The  accuracy  to  check 
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Newton-Raphson  convergence  was  0.1. 

All  these  computations  completed  one  iteration.  Further  iterations 
were  allowed  until  convergence  was  obtained. 

5.7   RESULTS 

Problem  A 

The  four  sets  of  initial  approximations  tried  in  this  case  are  listed 
in  Table  7.  Out  of  these  four,  three  sets  resulted  in  convergence  to  the 
solution  of  the  problem. 

The  convergence  rates  of  the  three  control  variables  and  the  six  state 
variables  for  problem  1A  are  shown  in  Figs.  15  through  23.   Fig.  24. shows 
the  optimal  profile  of  the  profit  function.   The  optimal  profiles  of 
Lagrange  multipliers  are  shown  in  Figs.  25  and  26.   The.  total  profit  was 
$79.75. 

In  order  to  visualize  the  convergence  rates  of  the  control  variables 
in  more  detail,  they  are  tabulated  in  Tables  13  through  15.         '  . 

Problem  3A  did  not  converge  to  the  optimal  solution.   It  encountered 
Newton-Raphson  convergence  problem  in  the  first  iteration. 

Some  modifications  were  made  to  study  the  convergence  problem. 

In  problem  4A,  T   was  changed  to  300,  with  other  values  remaining 
the  same.   This  problem  experienced  the  Newton-Raphson  convergence  problem 
in  iteration  2. 

Again  in  pr 

values  remaining  the  same,  the  advertisement  curve  became  nearly  discon- 
tinuous, but  it  did  converge  in  the  sense  that  there  was  little  difference 
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Fig.  15.   Convergence  Rate  of  Temperature  T  ,  Problen  1A. 
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Fig.  17.   Convergence  Rate  of  Advertisement  A,  Problem  1A. 
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Fig.  19.   Convergence  Rate  of  Concentration  y  ,  Problem  1A. 


81 


0.60 


0.55 


0.50 


£  0.A5 


0.40 


0.35 


-«._--_«_  1st  ITERATION 
.1.  6th,  optimal 


0.0      0.2      0.4      0.6 

TIME,  t 


1.0 


Fig.  20.   Convergence  Rate  of  Concentration  x  ,  Problem  1A. 
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Fig.  24.   Optimal  Total  Profit  Curve,  Problem  1A. 
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Table  13. 

Convergence 

rate  of  T  (t) 

in  problem  A. 

ite 

E  0 

1 

2      3 

4 

5 

6 

7 

8 

time 

0.0 

345.0 

366.3 

360.1  363.4 

361.9 

362.6 

362.3 

362.5 

362,4 

0.2 

345.0 

364.8 

355.8  360.8 

358.5 

359.6 

359.1 

359.3 

359.2 

0.4 

345.0 

364.4 

355.3  360.4 

357.9 

359.2 

358.6 

358.9 

358.7 

0.6 

345.0 

364.4 

356.2  360.5 

358.5 

359.5 

359.0 

359.2 

359.1 

0.8 

345.0 

363.7 

356.7  359.9 

358.4 

359.1 

358.8 

358.9 

358.9 

1.0 

345.0 

340.0 

340.1  339.9 

340.2 

340.2 

340.0 

340.0 

340.0 

Table  14. 

Convergence 

rate  of  T,(t) 

in  problem  A. 

iter  Q 

1 

2      3 

4 

5 

6 

7 

8 

time 

0.0 

345.0 

367.8 

364.8  366.0 

365.6 

365.7 

365.7 

365.7 

365.7 

0.2 

345.0 

365.8 

358.6  361.7 

360.4 

361.0 

360.7 

360.8 

360.8 

0.4 

345.0 

364.7 

356.2  360.3 

358.4 

359.3 

358.9 

359.1 

359.0 

0.6 

345.0 

364.3 

355.9  360.1 

358.0 

359.0 

358,6 

358.8 

358.7 

0.8 

345.0 

364.3 

356.7   360.3 

358.5 

359.4 

359.0 

359.2 

359.1 

1.0 

345.0 

340.0 

340.0  339.9 

340.2 

340.2 

340.0 

340.0 

340.0 

Table  15 

Convergence 

rate  of  A(t)  : 

Ln  problem  A. 

iter  Q 

1 

2      3 

4 

5 

6 

7 

8 

time 

0,0 

260.40 

358.80  366.40 

370.60 

369.90 

370.50 

370.30 

370.40 

0.2 

4.47 

4.65   4.68 

4.67 

4.67 

4.67 

4.67 

4.67 

0.4 

1.38 

1.70   1.72 

1.72 

1.72 

1.72 

1.72 

1.72 

0.6 

0.43 

0.67   0.69 

0.69 

0.69 

0.69 

0.69 

0.69 

0.8 

0.09 

0.22   0.23 

0.22 

0.22 

0.22 

0.22 

0.22 

1.0 

0.00 

0.00   0.00 

0.00 

0.00 

0.00 

0.00 

0.00 
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in  the  results  of  iterations  8  and  9. 

Problem  B 

All  the  three  initial  approximations  listed  in  Table  8  converged  to 
the  optimal  solution. 

Problem  3B  encountered  Newton-Raphson  convergence  difficulty  initially, 
but  using  the  same  grid  point  of  the  previous  iteration  as  the  starting 
value  in  the  Newton-Raphson  solution,  the  optimal  solution  was  obtained 
in  about  6  iterations.  The  optimal  profiles  of  the  six  state  variables 
and  the  three  control  variables  are  shown  in  Figs.  27,  28  and  28A.  The 
total  profit  in  this  case  was  $95.79. 

Problem  C 

Unfortunately,  for  both  the  initial  approximations  given  in  Table  9, 
the  problem  did  not  converge. 

In  problem  1C,  the  computer  experienced  exponential  overflow  while 
calculating  the  control  variables  in  the  first  iteration.  The  following 
values  of  the  state  variables  were  obtained  in  the  first  iteration 

1(1)  =  -12.62      S(l)  =  238.76 

All  other  values  were  reasonable. 

There  are  two  reasons  for  this  convergence  problem:   the  Newton- 
Raphson  convergence  difficulty  or  the  quasilinearization  difficulty.   In 
this  particular  case,  the  author  feels  that  it  was  the  Newton-Raphson 
convergence  difficulty. 
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In  problem  2C,  computation  was  not  possible  in  the  fifth  iteration 
because  of  the  Newton-Raphson  convergence  problem.   Results  of  the  fourth 
iteration  indicate  that  the  sales  goes  to  negative  in  the  initial  time 
and  then  rises  up  to  82.61  at  the  final  time.  Another  peculiarity  of 
this  problem  was  the  near  discontinuity  of  the  advertisement  curve.   A(0) 
is  -5639.0  and  A(0.01)  is  +90.54.   Such  a  sharp  change  of  the  variable 
with  time  can  make  both  the  Newton-Raphson  method  and  the  quasilinearizatlon 
method  useless. 

Examining  the  linearized  performance  Equations  (108)  and  (113) ,  we 


c.  can  make  the  problem  unstable.   In  this  case  c.  ■  0.0002.   The 
A  "  A 

author  feels  that  this  is  a  fairly  low  value  and  considering  that  Equations 


not  obtaining  a  solution  to  this  problem. 

In  spite  of  this  difficulty,  with  all  other  values  remaining  the  same, 
S(0)  was  changed  in  hope  of  obtaining  a  solution.  Initially  S(0)  was  0.1, 
and  the  following  values  of  S(0)  were  tried. 

a.  S(0)  =  15.0.   The  value  of  the  advertisement  was  negative, 
A(0)  =  -145.6.   This  case  encountered  the  Newton-Raphson  convergence 
difficulty  in  the  first  iteration. 

b.  S(0)  =  8.0.   The  value  of  the  advertisement  was  negative  and 
decreasing  very  rapidly.   A  solution  was  not  possible  because  of  the  Newton- 
Raphson  convergence  problem  in  the  second  iteration. 

c.  S(0)  =  4.0.  The  value  of  the  advertisement  was  negative  and 
decreasing  rapidly.   The  same  convergence  difficulty  was  encountered  but 
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this  time  in  the  fourth  iteration. 

It  can  be  seen  in  all  these  cases  that 

1.  The  value  of  the  advertisement  curve  was  nearly  discontinuous,  and 

2.  The  Newton-Raphson  method  caused  convergence  difficulty. 
This  leads  to  the  idea  of  increasing  c  . 

Problem  D 

c  =  0.01  and  S(0)  «  1.0 
a 

These  two  changes  were  made  in  the  parameters  of  problem  C.  Table  10 
shows  the  three  initial  approximations  used  in  this  problem.   Sets  2D  and 
3D  proved  to  be  good  guesses  and  convergence  for  these  sets  was  obtained 
in  about  5  iterations. 

This  problem  had  no  final  condition  on  the  inventory.  This  is  where 
it  differs  from  problem  A.   For  the  purpose  of  comparison,  the  detailed 
results  of  this  problem  are  given.   Figs.  29  through  37  show  the  conver- 
gence rates  of  the  three  control  variables  and  the  six  state  variables 
for  problem  2D.   The  convergence  rate  of  the  profit  function  is  given  in 
Fig.  38.   The  total  profit  was  $66.26.   The  advertisement  curve  as  can 
be  seen  in  Fig.  31  is  not  as  sharp  as  the  previous  ones.   It  is  a  mono- 
tomically  decreasing  curve. 

Set  ID  did  not  converge.   Examining  the  initial  approximations,  only 

^  n(t)  =  -0.5  could  be  a  wrong  guess,  but  the  final  solution  of  X,(t)  in 
b  ,U  o 

the  first  iterations  look  reasonable.   It  encountered  convergence  problem 
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Fig.  29.   Convergence  Rate  of  Temperature  T  ,  Problem  2C. 
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Fig.  31.   Convergence  Rate  of  Advertisement  A,  Problem  2D. 
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Fig.  33.   Convergence  Rate  of  Concentration  y  ,  Problem  2D. 
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Fig.  35.  Convergence  Rate  of  Concentration  y_,  Problem  2D. 
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As  a  slight  modification,  S(0)  ■  5.0  was  tried.  All  other  parameters 
remain  unchanged.   This  problem  encountered  the  Newton-Raphson  convergence 
difficulty  in  the  fourth  iteration.   Observing  the  optimal  profiles  of 
all  the  variables,  the  author  feels  that  this  was  the  most  stable  problem 
out  of  the  five  problems  tried. 

Problem  E 

Only  change  in  the  parameters  was  S(0)  ■  0.1,  the  other  parameters 
remain  unchanged.  The  values  of  the  initial  approximation  used  is  given 
in  Table  11. 

Set  IE  converged  to  the  optimal  solution.   The  optimal  profiles  of  the 
six  state  variables  and  the  three  control  variables  are  shown  in  Figs.  40, 
and  41.  The  total  profit  was  $65.98.   The  advertisement  profile  was 
very  sharp  again.   Obviously,  this  is  because  of  the  change  in  the  initial 
sales.  There  is  not  much  difference  from  the  other  profiles. 

In  order  to  get  an  overall  view,  A(0)  and  J  for  all  the  problems  solved 
are  compared  in  Tables  16  and  17. 

On  the  average,  this  problem  converged  in  6  iterations  with  3  digits 
accuracy.   For  9  iterations,  IBM  360/50  computer  took  about  16  minutes 
with  FORTRAN  IV  H  LEVEL  compiler.   The  computer  program  is  given  in 
Appendix  3. 

5.8   DISCUSSION 

The  results  for  all  the  problems  indicate  that  the  optimal  profiles 
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either  increasing  or  decreasing  slowly.   The.  optimal  profile  of  the  adver- 
tisement A  decreases  very  rapidly  making  it  almost  discontinuous.   Surpris- 
ingly, quasilinearization  did  not  encounter  any  trouble  with  this  type  of 
curve.   The  gradient  technique  [15]  and  the  second  variation  technique 
[13]  seemed  to  have  failed  because  of  this  curve. 

In  general,  convergence  was  obtained  in  5  to  6  iterations.   Tables 
16  and  17  give  the  comparison  of  the  convergence  rates  for  A(0)  and  J  for 
all  the  problems  solved. 

Comparing  the  optimal  curves  of  all  the  problems,  it  was  observed  that 
there  was  no  significant  difference  in  the  production  and  temperature 
profiles.   It:  may  be  concluded  that  a  change  in  certain  parameters  have 
little  effect  on  these  profiles.   However,  a  significant  change  in  the 
advertisement:  curve  was  noted.   This  can  be  explained  as  follows. 

With  certain  starting  values  of  production  and  inventory,  there  is 
a  definite  range  of  initial  sales  the  market  can  absorb  with  reasonable 
advertisement.   If  the  initial  sales  is  too  low,  the  market  needs  a  very 
high  advertisement  to  bring  up  the  sales  to  the  market  capacity.   On  the 
other  hand,  if  the  initial  sales  is  too  high,  the  market  needs  negative 
advertisement  to  bring  down  the  sales.  For  this  reason  initial  sales  was 
a  critical  value. 

Another  criti 

cost  of  advertisement  is  very  little.   From  the  cost  point  of  view,  heavy 

fluctuations  in  A  would  not  affect  the  optimal  solution  seriously.   Hence 

this  profile  was  observed  to  be  either  negative  or  discontinuous  or  unstable 

in  cases  where  c.  =  0.0002.   Stable  curves  were  obtained  with  c.  =  0.01. 
A  A 
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There  were  two  main  difficulties  in  solving  this  prohlem.   Convergence 
difficulties  in  1.  the  Newton-Raphson  method,  and  2.  the  Quasilinearization 
method.   It  was  because  of  the  former  difficulty  that  some  of  the  problems 
did  not  give  any  solution.  No  problem  failed  because  of  the  latter  diffi- 
culty.  If  a  better  method  for  solving  T  and  1^   could  be  used  the  author 
is  optimistic  that  all  the  problems  discussed  here  would  converge  to  the 
solution. 


CHAPTER  6 
CONCLUSION 

The  numerical  examples  presented  in  this  work  suggest  that  the  quasi- 
linearization  technique  may  be  a  useful  tool  for  obtaining  the  solutions 
of  nonlinear  mixed  boundary  value  problems.   Because  of  the  intimate 
association  between  the  boundary  value  problems  and  optimization  and 
control,  this  technique  is  also  a  useful  tool  for  solving  optimization 
problems  and  systems  analysis. 

The  object  of  this  work  has  been  to  illustrate  the  effectiveness  of 
this  method  in  overcoming  the  non-linearity  difficutly  in  two  point 
boundary  value  problems  encountered  in  optimization.  The  advantage  of 
this  approach  lies  in  its  rapid  rate  of  convergence,  provided  that  the 
initial  approximations  are  within  the  interval  of  convergence  of  the 
problem.   This  interval  is  fairly  large  for  a  number  of  problems.   Further- 
more, this  interval  can  be  enlarged  by  using  devices  such  as  data  pertur- 
bation. 

Convergence  rate  was  found  very  rapid  in  the  problems  solved.   This 
implies  computational  efficiency  in  terms  of  computer  time  for  a  prescribed 
accuracy.   This  is  an  important  advantage  of  this  method  over  other  optimi- 
zation techniques  such  as  the  gradient  techniques  where  the  convergence 
rate  is  slow  particularly  near  the  optimal  solution. 

The  convergence  of  this  method  is  contingent  on  the  choice  of  starting 
functions.   In  this  work  choosing  correct  initial  approximation  was  not 
difficult.   In  general  the  basic  physical  knowledge  of  the  system  is 
enough  to  make  a  correct  guess. 


114 


In  the  formulation  of  the  problem  the  control  variables  can  generally 
be  eliminated.   For  this  reason,  this  method  is  found  superior  to  other 
optimization  techniques  such  as  dynamic  programming  and  the  gradient 
techniques.  In  the  second  problem,  two  control  variables  could  not  be 
eliminated  from  the  performance  equations,  but  still  quasilinearization 
method  worked  well.  It  has  always  been  the  difficulty  in  solving  these 
control  variables,  that  some  of  the  problems  did  not  give  a  solution. 

This  method  was  observed  to  be  more  accurate  than  either  the  first 
variational  or  the  second  variational  techniques.   In  the  latter  methods, 
the  average  control  variable  is  used  throughout  the  calculation  of  one 
step  size.  For  a  fast  increasing  or  fast  decreasing  curve,  this  is  not 
likely  to  give  accurate  results.   Quasilinearization,  on  the  other  hand, 
uses  the  control  variable  for  the  calculation  at  the  same  point.  For 
this  reason  the  accuracy  is  higher  in  this  method. 

For  illustrative  purposes,  only  two  problems  have  been  considered. 
Obviously,  this  method  can  be  applied  to  a  variety  of  other  complex  problems 
arising  in  industrial  management  systems.  In  addition,  it  can  be  combined 
with  other  optimization  techniques  such  as  dynamic  programming  and  non- 
linear programming  to  optimize  various  topologically  complex  processes 
encountered  in  the  industry.  Recently,  this  method  has  been  proved  to 
be  an  efficient  tool  for  reducing  the  dimensionality  difficulty  in  dynamic 
programming . 
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NEWTON-RAPHSON  METHOD  OF  ROOT 

FINDING 

Basically,  this  method  is  the  Taylor  series 

expansion 

with 

second 

and  higher  order 

terms  neglected.   Expanding  f(u  , , )  around  u 

f<W 

-  f(u  )  +  (u  ..  -  U„)  f'(u  )  + 
n      n+1    n      n 

f(u  )  +  (u    -  uj  f'(u  )  - 

f(u  ) 
,1             n 

0 

n+1  "  n   f ' (u  ) 
n 

This  is  the  Netton-Raphson  equation  of  root  find: 

ing. 

Since  it  is 

essential  to  solve  two  equations  (for  T..  ; 

md  T 

)  simul- 

taneously  in  chaj 

iter  5,  we  derive  Newton-Raphson 

equations 

for 

such  case. 

f<\ 

TfVfd^,  T^)  +  (Tf1-^) 

+  (if1  -  T")  M  .  o 

3T2 

lit 

IT, 

SOf1, 

,  if1)  =  BdJ,  1$  +  (if1  -  T? 

+  (if1  -  A  2*  =  o 

L   3T2 

3T 

Rearranging  the  terms, 
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Tn+1  Jf  +     n+1  Jf  .  Tn  jf  +     n  Jf  _  £ (Tn     ^ 

1    at        2    3t2       hr1        3i2 


Since  t"  and  t"  are  known  (initial  approximation  needed  for  first 
iteration) ,  only  unknowns  in  these  equations  are  T  ,  and  T^  .  This 
iterative  procedure  is  carried  out  until  desired  accuracy 


if1  -  T*  <  E 


if1  -  T^<  e 


is  obtained. 

A  correct  choice  of  initial  approximation  should  be  emphasized.  Any 
error  in  this  selection  can  make  the  procedure  diverge. 
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APPENDIX  2 


JJOB  SHAH,RUN  =  CHECK,TIME  =  9,PAGES=100,LINES=5'j  119 

C 

C      THIS  PROGRAM  SOLVES  A  SET  CF  FOUR  DIFFERENTIAL 

C      EQUATIONS  TWO  POINT  SPLIT  BOUNDARY  VALUE  TYPE  USING 

C      SUPERPOSITION  PRINCIPLE  AND  RUNGE  KUTTA  TECHNIQUE 

C 

C 

c 

C      MAIM   PROGRAM 

C 

COMMCN  SX2,SX4,P,DT,A,B,CC,CA,AN,CI,AIM,C,J,Xl,X2,X3,X<i 
DIMENSION  X1(105),X2<  105),X3(105),XM105),SX1(1<3,102), 
lSX2(]-3,i02),SX3(l9,l02),SXM19,102),PXl!102),PX2(102), 
2PX3(102),PXM 102),H1X1(102),H1X2( 102 ),H1X3 1102 ),H 1X4 (102), 
3H2X1(102),H2X2(102),H2X3(102),H2X4!]02),PR1(105),PR2(105), 
4PR3C105)iPRFTU05)iA0VT1105) 

C 

C      READING  IN  DATA 

C 

100  F0RMAT(8F9.'t) 

READ  100,A,B,C,CC,CA,AN,CI,AIM 

101  F0RMAT(3F9.'O 
READ  101,DT,Wl,W2 

110  FORMAT! • 1VALUE  OF  THE  CONSTANTS') 

PRINT  110 
120  FORMAT!'  A= •  ,  F8 . 3 , ' 8= ' , F8 .3 , 'C= ' , F8. 3, 'CC= ' , F8.3, • CA= ' , 
lF7.3i,AN=,*F8.3,«CI='»F7.4,«AIM=«,F9.4,«Sf0CT)»«tF9.4i 
2«L6,Cm  =  '.F9.4.  'DT='  ,F6.3) 
PRINT  120,A,B,C,CC,CA,AN,CItAIH,Wl,W2.DT 
DO  130  1-1,101 
SX2«1,I)«W1 
SXM1.I )=W2 
130  CONTINUE 
140  FORMAT! •  XI ! 0 ) - ' , F10. A , '   X2  (0  )  =  '  ,  F  10.  4,  »   X3 ( 0) =« , F10.4, 

1«   X'.  (0)  =  *  ,  F10.4) 
HI  FORMATdH  ,I4,2X,4E20.7) 
C 

C      QUASILINEARIZATION  ITERATIONS  START 
C 

DO  3C0  J=2,19,l 
C 

C      PARTICULAR  SOLUTION 
C 

150    F0RMA1  UF10.41 

160  FORMAT! '-PARTICULAR  SOLUTION') 
PRINT  160 
P=l. 

RE  AD  1 50,  XI  (  1 )  ,  X2  (  1  >  ,  X3  ( 1  )  ,  X',  !  i  I 
PRINT  140,X1!1),X2! l),X3(  1) ,X4 !  1) 
CALL  RKT 
DC  170  1=1,101 
PXK  I)  =  X1  !  I  ) 
PX2I1  )  =  X2!  I) 
PX3!  I  )  =  X3!  I) 
PXM  I)=X4I  I) 


K 

170  CONTINUE 

120 

31 

C 
C 
C 

PRINT  141,  (  [,X1I  I),X2(l),X3(I),X4(I),I=l,101, 

20) 

HOMOGENEOUS  SOLUTION   FIRST  SET 

32 

P=0.C 

3  1 

180  FCr'.f'ATl '-HOMOGENEOUS  SOLUTION    FIRST  SET') 

34 

PRINT  180 

35 

READ  150,Xl(l),X2ll),X3(l),X4( 1) 

36 

PRINT  140,X1(1),X2(1),X3(1),X4(1) 

37 

CALL  P.KT 

3£ 

[10  ISO  1  =  1,101 

3<= 

H1XK  I  )  =  X1(  I  ) 

4C 

H1X2(1)=X2( 11 

41 

H1X3U )=X3(I  ) 

42 

H1X41 1 !=X4 ( 1  ) 

42 

190  CONTINUE 

44 

C 
C 
C 

PRINT  141,  (I, XII  I1.X2I I ),X3t I ) , X4 1  I ) , I  =  1 , 101 

,20) 

HOMOGENEOUS  SOLUTION    SECOND  SET 

45 

P  =  O.C 

4  6 

200  FORMAT ( '-HOMOGENEOUS  SOLUTION    SECOND  SET1) 

47 

PRINT  200 

4f_ 

READ  150, X 1 < 1 ) , X2 ( 1 ) , X3 ( 1 ) , X4 ( 1 ) 

4S 

PRINT  140, Xltll,X2ll),X3(l),X4(l  ) 

5C 

CALL  RKT 

51 

DO  210  1=1,101 

52 

H2X1II  I«X1(I) 

52 

H2X2( I)*X2(I ) 

54 

H2X3I  I1=X3(I  ) 

55 

H2X4I  I  )-X4(  I  ) 

5« 

210  CONTINUE 

57 

c 
c 
c 

PRINT  141, (I.XHI  ),X2(1 ),X3(I) ,X4(I1, 1=1,101 

,20) 

SOLUTION  OF  INTEGRATION  CONSTANTS 

5E 

220  F0RPATI2F9.4) 

5S 

READ  220.BBl.PB2 

6C 

B1=H61-PX3(101  ) 

61 

B2=EB2-PX4( 101 ) 

62 

CET  =  HX3(101)*H2X4( 101)-H1X4( 101)*H2X3! 101) 

62 

Al=tei*H2X4(101)-B2*H2X31101))/UET 

64 

A2=tP.2*HlX3(101)-Bl*lilX4(101)  1/DET 

65 

230  FORKATUH  ,  «  Al=  '  ,  F  9.4  , '     A2=  '.F9.4) 

66 

c 
c 
c 

PRINT  230.Al.A2 

RECOVERY  OF  SOLUTION       SUPERPOSITION  PRINCIPLE 

67 

DO  250  1=1,101 

It 

SX1  !  J,  I  )=PX1  (I  H  A1*H1XHI  )+A2*H2Xl<  I  ) 

6S 

SX2(J,T)=PX2(1 I+A1*H1X2«1 >+A2*H2X2( I ) 

7C 

SX3 i J,l )=PX3( I HA1*H1X3!I )+A2*H2X3(Il 

71 

SX4IJ.I  !  =  PX4U)  +  A1*H1X4(1)+A2*H2X4U) 

72 

250  CONTINUE 

26C  FORMAT!'  FINAL  SOLUTION       ITERATION  NO    '(Hi//l      121 
JJ=J-1 

PRIM  260,  J  J 
270  FORKATUH  ,I4,2X,'iE19.7,6X,E15.5,6X,E15.5) 
C 

C      CALCULATION  OF  CONTROL  VARIABLE  AND  PROFIT 
C 

P1R=C0 
P2R=C0 
P3R=C.O 

UO  2S0  1  =  1,  101 

ADVTI  I  )=(SXMJ, I )/(2.*CA) )* ( SX2 (J , I ) /AN- 1 . ) 
PRK I )=P1R+0T*(C*SX2< J, I ) ) 
PIR=PR1U) 

PR2 I  I )=P2R*DT*(CI*l lAIM-SXl(Jtl))**2)) 
P2R=PR2II) 

PR3I  I )-P3R:DT< (CA*SX2(J,I >*( ADVTI I )**2) ) 
B3ft=PR3l I ) 

PRFT<I)=PR1(I)-PR2U I-PR3I I ) 
*-iW-T-H-l  =  3.0 
260  CONTINUE 

PAINT  2  70,  ( t.SXK J.I ) ,SX2( J, I ) ,SX3( J, I ),SX<.< J, I) , ADVTI  I), 
1PRF-1  (I  1,1  =  1,101) 
290  FORMAT ('-TOTAL  PROFIT  ' , F15.3, 15X, 3F10. 3) 

PRINT  290,l>RFT<101),PRl(10l>,PR2(101>,PR3(101) 
C 

C      END  Cf  ONE  ITERATION 
C 

300  CONTINUE 
C 

C      END  CF  DO  LOOP  FOR  OUAS IL INEAR 1 1  AT  ION 
C 

STOP 
END 


5C0 


SUBROUTINE    KKT 

THIS  SUBROUTINE  I 
EQUATIONS    SIKULTA 

COMNCN  SX2,SX4,P, 
DIMENSION    XK105) 

1A2  (  1  C5)  ,A3(105)  ,A 

28M1C5),C1  (105)  ,C 

3D2(lC5),rj3ll05),D 
DO    5C0    1^1,  100 
«=SX2(J-1 ,1 ) 
N«SX4(J-1, I ) 
TAM--1 
T«TA*DT 

AH I >»DT*(P*A+P*B 
B1A*CC*V+W*IV**2> 

lH»lV**3l/t2.*CA*< 
B2A=CC+2.*V*W/(CA 

1/(2. *CA<  IAN**2) ) 
B3A*V**2/(CA*AN)- 
B1U  >«DT*(P*B1A-P 
CL! ! »=DT*(P*2.*C1 
DlA=C-W*CC+3* *(W* 

1-IW**2)*V7 ICA*AN1 
D2A=3.*V*(W**2)/( 
03A-V,/(2.*CA)-CCi 

UCA*AN) 

1«>2.*CC*V/AN 
01 ( I )=DT*tX3(  I  M 

1D3A) 
A2( I >«DT*tP*A+P*B 
B2(r  I=0T*(P*B1A-P 

UX'.l  I  !  +  D  1(11/2.  )* 
C2(I»*0T*(P*2.*CI 
D2( I )^ET*(  !X3(  ID 

1  +  B1! I  1/2. )*D2A+(X 
A3(I  >=DT*(P»A-*P*B 
B3U)«DT*(P*BIA~P 

UX41IDD2U  1/2.  )* 
C3(  1  »eOT*{P*2.*Cl 
D3(H=BT*I  (X31D  + 

1*B2(  D/2.DD2AHX 
A4(  I )-DT*(P*A+?*E 
BMI  )  =  DT*ir>*BlA-P 

1UM  ID03U  )  !*33A 
CI  I  )=DT*(P*2.*CI 
DMl  l*DT*<  tX3UD 

IB3!I i  >*D2A  +  (X4U ) 
Xl( I tl)=Xl  tIMIAl 

X2(  r»l)>x2(  [)+(Bl 

X3!I  +  1)  =  X3(D  +  (C1 
X'.l  I-H)  =  X4  I  I  )  +  (Dl 
CONTINUE 
RETURN 
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S    USED    TO    INTEGRATE    THE    FOUR    LINEARIZED 
NEOUSLY    BY    RUNGE    KUTTA    METHOD 

DT,A,,f.,CC,CA,AN,CI,AIM,C,  J,X1,X2,X3,X4 
,X2il05),X3il05),X<.(105),Al(105), 
M1051.B1  (105)  ,B2(  105  1,63  D05), 
2(105), C3(105),C4 (105), 01(1 05), 
M105DSX2I  19,102) ,SXM 19, 102) 


*T-X2(  ID 

/<CA*AN)-W*V/(2.*CA)-CC*(V**2)/AN- 
AN**21  ) 
*AN)-W/(2.*CA)-2.*CC*V/AN-3.*W*(V**2) 

V/(2.*CA)-V**3/ (2.*CA* (AN**2) ) 

*V*B2A-P*W*B3A+X2(  (DB2A  +  XM I )*B3A) 

*AIH-2.*CI*X1< I ) ) 

*2)*(V**2I/i4.*CA*(AN**2>)+W**2/(4.*Cit) 

♦Z.*CC*V*h/AN 

2.*CA*(AN**21  )-W**2/(CA*AN)+2.*CC*U/AN 

}.*W*(V**2 ) / (2.*CA* ( AN**2 ) l-2.*W*V/ 


P*V*D2A-P*W*D3A*X2tt>*D2A+X4(I>* 


ID 

V*E 
3A) 
A  I  <■; 
Ml 
(!) 

it* 
\r*B 

3A) 

tin 

21  i 

in 

(D 

\t*B 


DT/2.)-(X2(I DBK I)/2.) ) 
?.A-P*W*P.3A+(X2(  I)+B1(I)/2.)*B2A< 

-2.*CI*(X1  U  DAltl)/?.)) 

1/2.  )+P?DLA-P*V*D2A-P*l-;'iD3A+(X2(  I  ) 

i01 (I)/2.)*03A) 

Dl/2. !~!,VtIDB2(U/2.D 

2A-P*W*B3*+lX2(  I  )+B2(I  1/2.DB2A  + 

-2.*CI*!XMIDA2(  l)/2.  D 

1/2.  l*P*DlA-P*V*D2A-f»«H*D3A+tX2lIJ 

+  D2! I )/2. !*D3A) 

DT)-CX2<  I)+B3(  ID) 

2 A -P*H.*B3A+IX2(  I  )+83(  I)  )*B2A  + 


IK-2.*CI*(X1(I 1+A3ID ) ) 
:3(  I! D-P*D1A-P*V*02A-P*W*D3A*  !X2(I 1  + 
»D3( II DD3A1 

!  I  }*2.*A2(I)«2.*A3(X)-»A4(I)  )/6. 
! I l+2.*B2(I)+2.*B3(1>+B4(I) 1/6. 
(I)+2.*C2U)+2.*C3(I>+C4U>)/6. 
! I D2.*D2( I l+2.*D3l I) »04( I  1 1/6. 
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APPENDIX  3 


$Jt 

c 

c 

c 

c 

c 

c 

c 


SHAH,RUN=CHECK,TIME=9,PAGES=100,LINES=55 
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THIS  PROGRAM  SOLVES  A  SET  OF  TWELVE  DIFF.  EQUATIONS  TWO  POINT 
SPLIT  TYPE  USING  THE  SUPERPOSITION  PRINCIPLE  AND  RUNGA_KU1 TA 
TECHNIQUE. 


MAIN  PROGRAM 

COMMCN  Xl,YltX2,Y2,AItQ,Zl,Z2,Z3,Z4,Z5 

lGA,Ce,FA,R,AQ,Vl,V2,Y0,EB,ANtCA,CItAIM 
2CT,P,W1,W2,K3,W4,TEM1,TEM2, J,SAQ,SZ6 

DIMENSION  X1(105),Y1(105),X2(105),Y2(1 
1Q(105),Z1(105),Z2(105),Z3(105),Z4(105) 
2Z6(1C5),SX1( 10,102),SY1UO,102),SX2(10 
3SAI(10,102),SAQ(10,102),SZ1(10, 1021.SZ 
4SZ3(10,102),SZ4< 10,102),SZ5I10,102),SZ 
5PX1(102),PYK102),PX2(102),PY2(102) ,PI 
6PZ1(102),PZ2(102),PZ3I102),PZ4(102),PZ 
7X0(1C5),Y0(105),TEM1(10,102),TEM2(10, 1 

DIMENSION  H1XK  102)  ,111  Yl(  102),  HI  X2(  102 
1H1AI(102),H1Q(102),H1ZK102),H1Z2(  102) 
2H1Z4(102),H1Z5(102),H1Z6(102),H2X1(  102 
3H2X2(102),H2Y2(102),H2AI(102),H2Q(102) 
4H2Z2U02),H2Z3(102),H2Z4(102),H2Z5(102 
5H3X1(102),H3Y1(102),H3X2( 1 02 ) , H3Y2 ( 102 
6H3Q(102),H3Z1(102),H3Z2(102),H3Z3(102) 
7H3Z5I 102),H3Z6( 102) ,H4X1 ( 102 > ,H4 Yl ( 102 
8H4Y2(102),H4AI ( 102 ) , H4Q ( 102 ) , H4Z 1 1  102) 
9H4Z3U02)  ,  H4Z4U02)  ,H4Z5  (  102  )  ,  H4  Z6  (  102 
1H5YK102),H5X2(102)  ,H5Y2  (  102  )  ,H5AI  ( 102 
2H5Z1(102),H5Z2(  102 ) ,H5Z3 ( 102) ,H5Z4 ( 102 
3H5Z6I102) ,H6X 1(102) , H6Y1 ( 102 ) , H6X2 ( 102 
4H6AII102),H6Q(102),H6Z1(102),H6Z2(102) 
5H6Z4I  102),H6Z5(  102)  ,H6Z6(  102) 

DIMENSION  ADV(105),PR1(105),PR2( 105I.P 
1PR5U05),PR6(105),PRFT(105) 

READING  IN  DATA 


,Z6,TP,TQ,DT,X0, 
,CB,CC,CO,C,AC, 


05), All 

,Z5(105 

,102) ,S 

2( 10, 10 

61 10,10 

(102)  ,P 

51 102), 

02) 

)  ,H1Y2( 

,H1Z3(1 

),H2Y1( 

,H2Z1(1 

) ,H2Z6( 

),H3AI( 

,H3Z4(1 

)  ,H4X2( 

,H4Z2( 1 

),H5X1  ( 

),H5QI1 

) ,H5Z5( 

1.H6Y2I 

,H6Z3(1 


105), 

), 

Y2(10,102), 

2), 

2), 

0(102) , 

PZ6(  102)  , 

102), 

02)  , 

102), 

02), 

102), 

102)  , 

02), 

102), 

02)  , 

102), 

02)  , 

102), 

102), 

02), 


R3( 105), PR4! 105), 


401  FORMAT(2E10.3,2F7.1,2F5.1,5F4.1,F5.1,F3.1 ) 
READ  401,GA,CB,EA,EB,TP,T0,R,AQ,V1,V2,AIM,AN,C 

402  F0RMAT(F7.5,F4.2,F7.5,5F5.3,4F6.3,F5.1) 
READ  402,CT,DT,CA,AC,CR,CC,C0,CI , Wl , W2 ,W3 , W4 ,T IM 

406  F0RKAT11H1, 'VALUE  OF  THE  CONSTANTS'! 
PRINT  406 

403  F0SMAT(4H-GA=,E10.3,'   GB=',E10.3,'   EA=«,F7.l,' 
1F7.1,'   R=',F':.l,'   FLOW  RATE=  '  ,  F5.  1  ,  '   V1  =  ',F4.1 
2'   V2=',F4.1,'   MEAN  INV.= ' , F4. 1 , •   MEAN  TEMP= • , F6. 2 ) 

PRINT  403,GA,GB,EA,EB,R,AQ,V1,V2,AIM,TIM 

404  FORMAT!  1H-,  »N=', F5.lt'   C-NF3.lt'   CT=',F7.5,< 
1F4.2,'   CA=',F7.5,'   C1=',F5.3,'   C2=',F5.3,'   C2 
2F5.3,'   CQ=',F5.3,'   CI=',F5.3,'  XO ( T ) = ■ , F6. 3 , 

3'  YO(T)=' ,F6.3) 
PRINT  404, AN,C,CT,DT,CA,AC,CB,CC,C0,CI,K1,W2 

405  FORMATUH-, 'T1=',E12.4,  ■   T2=',E12.4,'   00  (  T  )  =  •  ,  E  10.  2  , 


1  i 


DT='  , 

:3=>, 


1'       L6(T)=' ,E10.2)  125 

PRINT    405,  TP,TQ,W3,W4 
129    FORMATI1H    ,  14,  2X,  12E10.2  ) 
110    F0Rt:AT(lH-,'Xl(0)  =  ',F4.2,'Yl<0)  =  ',F4.2,,X2(0>  =  'fF4.2, 

l'Y2(0)=',F4.2,'I  (0)  =  '  ,  F4.  2  ,  ■  CM  0  )  =  •  ,  F4.2 , ' LI (  0  )  =  '  ,F4.2, 

2'L2(0)=',F4.2,'L3(0)=,,F4.2,'L4(0)=',F4.2,'L5(0)=',F4.2, 

3'L6(0)=',F4.2> 

725  F0RMATI4E15.5) 

726  FORMAT  UH  ,  'Tl  =  ' ,  El 5.  5,  'T2=  •  ,E15.5,  '0=  ■  .E15.5,  'L6='  ,E15.  5) 
DO  162  1=1,101 

TEMU1,I)  =  TP 

TFM2I 1,T)=TQ 

SAC! 1 ,! >=U3 

SZ6I1.I )»W* 
162  CONTINUE 
C 

C      QUASH  I  NEAR  RATION  ITERATIONS  START 
C 

00  3C0  J  =  2, 10, I 
C 

C      PARTICULAR  SOLUTION 
C 

601  FORMAT ( 12 (F5.2) ) 

600  F0RMATI1H-, 'PARTICULAR  SOLUTION') 

PRINT  600 

P=L. 

READ  601,X1(1),Y1(1),X2(1),Y2(1),AI ( 1 ) ,Q ( 1 ) , Z 1 ( 1 ) , 
lZ2(l),Z3(l),Z4(l),Z5(i),Z6(l) 

PRINT  110,X1(1),Y1(1),X2(1),Y2( t ),AI(1),Q(1  ),Z1(  1), 
2Z2(l).Z3(l ),Z4(1),Z5( 1),Z6( 1) 

CALL  RKT 

DO  2C1  1=1,101 

PXK  I  )  =  X1(  I  ) 

PYK  l)  =  Yl(  I) 

PX2( I)=X2(I ) 

PY2I I )  =  Y2( I ) 

PI ( I )  =  AI ( I  ) 

PQ(I )=QII) 

PZUI)  =  Z1(I) 

PZ21 I)=Z2(I) 

PZ3! I)=Z3( I) 

PZ4( I)=Z4( I) 

PZ5I I )=Z5(I ) 

PZ6(I)=Z6(I) 
201  CONTINUE 

PRINT  129,  I  I, XII  I ),Y1(I),X2( I),Y2( I ),AI( I),Q(I),Z1( I), 
1Z2(I),Z3(I),Z4(II,Z5(I),Z6(1),I=1,101,20) 
C 

C      HOMOGENEOUS  SOLUTION   FIRST  SET 
C 

P=O.C 

602  FORMATUM-'HOMOGENEOUS    SOLUTION    FIRST    SET') 
PRINT    602 

READ    601,X1(1),Y1 (1),X2(1) ,Y2( 1) ,AI ( I) ,Q(1),Z1 (1) , 
1Z2(1),Z3(1),?4(1),Z5(1),Z6(1) 


PRINT  llO,Xlll),Yl(l),X2(l),Y2d),AI(l),Q(l),Zltl), 
2Z2(1).Z3(1),Z'<(1),Z5(1),Z6(1) 
CALL  RKT 
DO  2C2  1=1,101 
H1X1(I)=X1(I ) 
H1Y1(I)=Y1 (I) 
H1X2(I)=X2(I ) 
H1Y2(I)=Y2(I) 
H1A1  (I)=AI  !I) 
H1Q1I)=C( I ) 
HIZ1(I)=Z1(I) 
H1Z2(I)=Z2(I) 
HlZ3Ul=Z3tI) 
HIZMI)=Z4(I) 
H 1 7.  C.  (  I  )  =  Z  5  ( 1  ) 
H1Z61  I  )  =  Z6II  ) 
202  CONTINUE 

PRINT  129, (I,X1  (I).Yl (I),X2(I ),Y2U l.AIt ! ),Q( I ) , Z  1  ( I), 
1Z2  (  I  ),  Z3(  I  ),ZM  I  ),Z5(  I  ),Z6(  11,1  =  1.101,20) 

HOMOGENEOUS  SOLUTION  SECOND  SET 

P  =  0-0 
603  FORMATdH-, 'HOMOGENEOUS  SOLUTION  SCCOND  SET<) 

PRINT  603 

READ  601,X1(1),Y1(1),X2(1),Y2(1),A1 (1) ,0(11,71(1), 
172(1),  Z3(1),7M1),Z5(1),76(  I) 

PRINT  110,Xl(l),Yl(l),X2(l),Y2d),AId),Qd),Zl(l), 
2Z2(  l),Z3<l),Z4(l),Z5(  ll.ZMl) 

CALL  RKT 

DO  2C3  1=1,101 

H2XK  I)  =  X1(I  ) 

H2Y1(I)  =  Y1(  I) 

H2X2(I)=X2U) 

H2Y2(I)=Y2(I) 

H2AH I)=AI ( I ) 

H2Q(I)  =  Q(I  ) 

H2ZKI)=Z1(I) 

H272(I)=Z2(I> 

H2Z3(I)=Z3(I) 

H2ZMI)  =  Z<i(I) 

H2Z5tl)=Z5d) 

H2Z6(  U  =  76(  I  ) 
203  CONTINUE 

PRINT  l?9,(I,XKI>,Yld),X2(I),Y2(l),AI(l),Q(I),Zl(T), 
lZ2(I),Z3(I),Z4!I),Z5(I),Z6(t),I=l,101,20) 
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HOMOGENEOUS  SOLUTION  THIRD  SET 

P  =  CC 
604  FORMATdH-, 'HOMOGENEOUS  SOLUTION  THIRD  SET') 
PRINT  604 

READ  601,X1(1),Y1(1),X2(1),Y2( I) ,AI ( 1 ) ,Q ( 1 ) , Z 1 (1 ) , 
1Z2(1),Z3I  1),Z'1(1),Z5(1),Z6(1) 
PRINT  110,Xld),Yl (1),X2(1),Y2(1),AI(1),C(1),Z1(1), 
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2Z2(l).Z3(l),Z4(l),Z5(l),Z6(l) 

CALL  RKT 
00  2C4  1=1,101 
H3XK  I)  =  X1(  I  ) 
H3Y1 ( I)=Y1(I) 
H3X21I)=X2(I ) 
H3Y2I I)=Y2(I  ) 
H3A11  I)  =  .M  (  I) 
H3QII)=Q(I  ) 
H3Z1(I)=Z1(I) 
H3Z2I  I  )  =  Z2(I ) 
H3Z3I [)=Z3(1 ) 
H3Z4(I)»Z4(tJ 

H3Z5(  n  =  z5(i ) 

H3Z6I l)=Z6(I ) 
204  CONTINUE 

PRINT  129, (I, XI (I ),Y1(I),X2( I ) , Y2 ( I ) , AI ( I ) , Q ( I ),Z1(I), 
1Z2(I),Z3(I),Z4(I),Z5(1),76(I),I=1,101,20) 
C 

C      HOMOGENEOUS  SOLUTION  FOURTH  SET 
C 

P=0.0 
605  FORMATIIH-, 'HOMOGENEOUS  SOLUTION  FOURTH  SET1) 
PRINT  605 

READ  601.XK  1  )  ,  Yl  ( 1 )  ,  X2  (  1 )  ,  Y2  (  1 )  ,  Aid),  0(1)  ,Z1(L)  , 
1Z2(1),Z3(1),Z4(1),Z511 ),Z6(l> 

PRINT  110, X1(1),Y1(1),X2(1) ,Y2(1),AI(1),Q(1),Z1(1 ), 
2Z2(U,Z3(1),Z4(1),Z5(  l),Z6(  1) 
CALL  RKT 
DO  2C5  1=1,101 
H4X11I  )  =  X1(I  ) 
H4YK  I)  =  Y1(I  ) 
H4X2II)=X2(I ) 
H4Y2U )=Y2(I ) 
H4A I ( I )  =  A  I ( I ) 
H4Q(1)=0(I) 
H4Z1(I)=Z1(I) 
H4Z2I  I)  =  Z2(I) 
H4Z3(I)=Z3(I ) 
H4Z4(I)=Z4(I) 
H4Z5(I)=Z5(I) 
H4Z6I I)=Z6(I) 
205  CONTINUE 

PRINT  129, (I, XI (I) ,Y1(I),X2(I),Y2(I),AI(I),0(I),Z1(I), 
1Z2(I),Z3(I),Z4(I),Z5( I ) ,Z6( I ) , 1=1 , 101 , 20) 
C 

C      HONCCENEOUS  SOLUTION  FIFTH  SET 
C 

P=O.C 
606  FOR.vATdH-, 'HOMOGENEOUS  SOLUTION  FIFTH  SET') 
PRINT  606 

READ  601,X1(1),Y1(1),X2(1),Y2(1),AI(1),Q(1),Z1(1), 
1Z2(1),Z3(1),Z4(1),Z5(1),Z6(1  ) 

PRINT  110,X1(1),Y1(1),X2(1),Y2(1), A  I ( 1 ) , 0 ( 1 ) , Zl ( 1 ) , 
2Z2(1),Z3(1),Z4(1),Z5(1),Z6(1) 


CALL.  RKT 

00  2C6  1*1,101 

H5XMI)  =  X1(I) 

H5Y1(I)=Y1(I) 

H5X21 I)=X2(I  ) 

H5Y2( I)=Y2  1  I  ) 

H5AU  I  )  =  A  I  (  I  ) 

H5QtI)=Q( I  ) 

H5ZUI)=Z1(I) 

H5Z2(I)  =  Z2(  I  ) 

H5Z3(I)=Z3(I) 

H5ZM I)  =  ZMI) 

H5Z51I J=Z5(I) 

H5Z6I I)=Z6(I) 

206  CONTINUE 

PRINT  129, (I,X1(I),Y1(I),X2(I),Y2(I),AI(I ),Q(I),Z1(I) 

c 

c 
c 

U2 1 1  ),Z3(  I  ),Z4(  I), Z5U),Z6(  I  1,1  =  1,101,201 

HOMOGENEOUS  SOLUTION   SIXTH  SET 
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P=0.0 
607  F0RFATI1H-, •HOMOGENEOUS  SOLUTION  SIXTH  SET') 

PRINT  607 

READ  601,X1<1),Y1(1),X2(1).Y2(1),AH1),Q!1>,Z1(1  ), 
172  I  1 )  ,  Z3  <  1  )  ,  Z4  ( 1 )  ,  Z5  I  1 )  ,  Z6  I  I  ) 

PRINT  110,X1(1),Y1(1),X2(1),Y2I1),AI(1),0(1),Z1(1)I 
2Z2(1)IZ3(1),Z4(1),Z5(1),Z6(1) 

CALL  RKT 

DO  2C7  I--1.101 

H6X1!  !  )  =  X1  (I  ) 

H6Y1!I)=Y1 (  I  ) 

H6X2(I)=X2(1  ) 

H6Y2( I)=Y2(I ) 

H6AI(I)=AI ( I) 

H6Q( I)  =  0( I ) 

H6Z1(I)=Z1(I) 

H6Z2(I)=Z2(I) 

H6Z3U)  =  Z3(I) 

H6ZMI)  =  ZMI  ) 

H6Z5(I)=Z5(I ) 

H6Z6I I )=Z6(I ) 
2C7  CONTINUE 

PRINT  129, (I.XK I),YHI),X2(I),Y2(I),AI(I),0(I),Z1(I), 
1  Z2(  I  ),Z3(I),Z4 1 1 ),Z5( I), Z6< 11,1  =  1,101,20) 

SOLUTION  INTEGRATION  CONSTANTS 

DIMENSION  E(6),A(36),BB(6) 
N=6 
150  F0RHAT(6fF5.21 ) 

READ  150, (BC( I ), 1=1,6) 
800  FORMTIIH  , 'FINAL  CONDIT IONS • , 10X, 6F15.4 > 
PRINT  800, (BB1I) , 1=1,6) 
B(l)=BB!l)-PZ5(101) 
B(2)=B8(2i-PZl(101) 


B(3)=8B(3)-P7.2(101)  129 

B(4)=BB(4)-PZ3(101) 

B(5)=BB(5)~PZ4(101) 

P(6}=BB(6)-PZ6(10l) 

A( 1)=H1Z5(101) 

A(2)=H1Z1(101) 

A(3J=HIZ2( 101) 

A(4)=H1Z3<  101) 

A(5)=H1Z4(  101) 

A(6)=H1Z6( 101) 

Af7)=H2l5i 101) 

A(8)=H2Z1(101) 

A(9)=H2Z2( 101) 

A(1C)=H2Z3(101) 

A{  11)=H2ZM101) 

A(12)=H2Z6( 101) 

A(13)=H3Z5(101) 

Al 14)=H3Z1(101) 

A( 15)=H3Z2I101) 

A( 16)=H3Z3(101) 

A(17)=H3Z4  (  101) 

A(18)=H3Z6(101) 

A( 19)=H4Z5(101) 

A(20)=H4Z1<101) 

A(21)  =  H',Z2(  101) 

A(22)=H4Z3(101  ) 

A(23)=H4Z4  (101) 

A(2M-H<tZ6(  101) 

A(25)=H5Z5(  101) 

A(26)=H5Z1U01) 

A(?7)=H5Z2(101> 

A(28)=H5Z3(101) 

A(29)  =  H5ZM101) 

A(30)=H5Z6(101) 

AC31)=H6Z5( 101) 

A(32)=H6Z1 (101) 

A(33)=H6Z2(101) 

A(34)=I16Z3(  101) 

A(35)  =  H6Z'.(101) 

A(36)=H6Z6(101) 

KS=0 

CALL    SIKQ( A,B,N,KS) 
151    F0RPAT(1H-,'A1=',F15.5,,A2=',F15.5,,A3=,,F15.5,  •  A4=  •  , 
1F15-5,,A5=,,F15.5, •A6=',F15.5) 

PRINT    151, (B(l 1,1=1,6) 
C 

C  RECOVERY    OF    SOLUTION  SUPERPOSITION    PRINCIPLE 

C 

00    160    1=1,101 

SXIU,I)  =  PXI(IMB(1)*H1X1(I)+B(2)*H2XI(I)+B(3)*H3X1(  I)  + 
ie(4)*H4XH I )+B(5)*H5XK I )+B I6)*H6Xl ( I ) 

S¥i(  J,I)=PY1(I  )+B(  1KH1YK  I  )+B(2!*H2Yl(I  )+B(3)*H3Yl(  I  )•» 
18(4)*H<m(  I  )+B(5)*H5Yl  (I  )+B(6)*H6Yl  (  I  1 

SX2( J,I)=PX2(1 )+B(l)*HlX2(I )+B(2)*H2X2( I ) +B ( 3 ) +H3X2 ( I )+ 
1B(4)*H4X2U  )  +  B(5)*H!>X2(I)  +  B(6)*H6X2(I  ) 


SY2(J,I)=PY2(I)+B(1)*HIY2<I  ) +B  (  2  )  *l  !2Y?  (  I  ) +  B ( 3) *H3Y2 < I )  I 
1BI4)*H4Y? ( I  )+B(5)*H5Y2(l  )+B(6)*H6Y2II) 

SAI(J,I)  =  PI<I)  +  B(l)*lilAI  (I)+B<2)*H2AI ( I ) +BI 3 )*H3AI  (  I H 
1B(4)*H<.AI  I  I  )+B(5)*H5Al  (  I  1+B(6)*H6AI  (I  ) 

$AQ(J,n---PQ(I>  +  H(l)*HlQ(IHB12)*H2Q(I)+Bl3)*H3Q<I)  + 
1B(4)*H4  0(I  )4B15)*H5Q(I)+Btt)*H60(  I  ) 

SZ1(J,T)  =  PZ1U  HB(  11*11171(1  )+B<2)*H2ZKI)+B(3)*H3ZlU)  + 
IB(4)*H421  (  I  )+BI5T*H5Zl(l  )+B(6)*H671(  I  ) 

SZ2<  J,  I  l*.pZ2m+BU)*HXZ2U  )  +  B  12  )*II2Z2  ( I  ) +B  (  3  )*H3Z2  <  I  »♦ 
1B(4)*H4Z2I I )+B<5)*H5Z2( I ) +B(6)*H6Z21 1  ) 

SZ3(J.1)=PZ3(I)+B(1)*H1Z3(I ) +B ( 2 ) *H2Z3 (l) +B ( 3 )*H3Z3 I  I ) i 
1B(4 )*H4Z3(I  )  +  B(5)*H5Z3(  I  ) +B < 6 ) *H6Z3 ( I  ) 

SZ4U,I)=PZ4(1)+B(1)*H1Z4(I)+B(2)*H2Z4(I  ) +B ( 3 ) *H3Z4 (I ) + 
ie!4)*H4Z4I  I  )+Bl5)*H5Z4(  I  )+B(6)*H6Z4(l  ) 

S/.5<  J,I)  =  PZ5(I  )+B(l)*HlZMI  > +B < 2 ) *H2 Z5 ( I  HB(3)*H3Z5(  1)  + 
1B(4)*H4Z5(  I  )+B(5)*H5Z5(  I  ) +  B  (  6  )  *H6Zrj  (I  ) 

SZ6(J,I)=PZ6<I)+B(1)*H1Z6(I)+B(2)*H2Z6(I)+B[3)*H3Z6U)+ 

1B(4J*H4Z6( I )+B15)*H5Z6(I ) +B (6 )*H6Z6(  I  ) 
160  CONTINUE 
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PRINTING  THE  FINAL  SOLUTION 


161 


ITERATION  NO  • ,13) 
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F0RMTI1H1,  'FINAL  SOLUTION 

JJ=J-1 

PRINT  161, JJ 

FORMATUH  ,  I'.,3X,6C18.5) 

PP.  I  NT  131,  !  I.SX1I  J,  D.SYIIJ,  I)  tSX2l  J,I  ),SY2(J,  I), 
1SAIIJ, Il.SAQU, I), 1=1,101 ) 
132  F0RMATUH-|« ADDITIONAL  STATE  VARIABLES*,//) 

PRINT  132 

PRINT  131,  (I,SZ1IJ,I),SZ2(J,  D.SZil  JtI)tSZ4  !  J,  I), 
1S75(J,I),S76(J,1 ), 1  =  1, 101) 

CALCULATION  OF  CONTROL  VARIABLES  AND  TOTAL  PROFIT 


4,17X,«Tl<=,,E14,6t3X,«T2»«, 

5) 


164  FGRMATIIH  ,  14  ,  5X  , '  ADVT=  •  ,  E  12 
1E14.6.12X, 'TOTAL  PROFIT', F10 
P1R=G.0 
P2R=0.0 
P3R=C.O 
P4R=C.O 
P5R=C.O 
P6R=C.O 

DO  5C0  1=1,101 

AOV( l)=(SZ6(J,I)/(2.*CA) )*( ( l./AN)-Il./SAQ(J,l>) ) 
PR1(1)=P1R+AC*CQ*SAQ< J, I l*DT 
PlR=PRli I ) 

PR2C1 )=P2R+CB*AQ*SX2( J,I)*DT 
P2R=PR2 I  I ) 

PR3I I)=P3R+CC*AQ*(1.-SX2( J, D-SY2I J,l ) )*DT 
P3R=PR3( I! 

PR4(  I)  =  P'iR  +  CI*(  (AIH-SAK  J,  I)  )**2)*DT 
P4R=PR4( I) 
PR5(I)=P5R+CA*(ADV(I)**2)*(SAQ(J,I)**2)*DT 

P5R=PRS! I ) 


131 
PR6I !)=P6R+CT*( (1 IK-TEK2  1 J- 1,1 ) )**2+(TEMl IJ-1,1) 

1-TEM2U-1.  r  !  1**2)*DT 

P6R=:pR6( I ) 

PRFTI  I)  =  PR1(  I1+PR2I  IHI'K.M  1  )-PR4(  I  )-PR5(  I  I-PR6I  I  ) 

500  CONTINUE 

502  FORMAT (1H-, 'TOTAL  PROFIT'' ,F15.7, 10X,6F15.5) 

PRINT  502,  PRFTI 101  ) , PR1 ( 101 ) , PR2 I  101 ) , PR3 (  101)  ,PR1(10l) 
I,PR5( 101),PR6(101I 

163  FORKATUH1, 'VALUES  OF  THE  CONTROL  VARIABLES') 
PRINT  163 

501  FORMATUH  ,'    ADVT  FOR  PREVIOUS  ITERATION         TEMP 

1  FOR  CURRENT  ITERATION' i 12X, ' PROFIT  FOR  CURRENT  ITERATION') 

PRINT  501 
C 

C      CALCULATION  OF  TEMP  1  AND  TEMP  2  BY  NEWTON  RAPHSON  METHOD 
C 

DIMENSION  TMK202)  ,TM2(202) 

IF  (J.GE.3)  CO  TO  905 

TKK1I--350.0 

TM2( 11=350.0 

GO    TC    906 
903    IH1I1 )=TEM1( J-l, 1 ) 

TM2(1)=TEM?(J-1, 1) 
906    CONTINUE 

DO    165    1=1,101 

DO    166    N=1.201 

D1T1A=USZ1(J,IJ-SZ2{J,1))*6A*SX1IJ,I)*EA)/(R*(TMHNI**4)I 

D1T1D=(EA/R-2.*TM1(N) )*EXP (-EA/ I R*TM1 I N ) ! ) 

DlTlf>tSZ2(J,I  )*Gt>*SYlt  J,  I  I*E81/(R*(TM1  IN)**'. )  ) 

DltlC=(EB/R-2.*TMl(NI )*EXPt-EB/(R*TMl IN)  !  ) 

D1T1=01TIA*01T16+01T1C*D1T10-4.*CT 

D1T2=2.*CT 

D2T1=2.*CT 

D2T2A=( ISZ3IJ, I 1-SZ4  1 J, I ) )*GA*SX21J( I ) *EA ) / I R* 1 1 M2 IN) v< 4 )  ) 

D2T2P=(EA/R-2.*TM2(N) )*EXP (-EA/ t R*TM2 ( N ) ) ) 

D2T2C= ( SZA (J , I ) *GB*SY2 ( J, I J*EB)/ (R* t TM2IN ) **4 ) ) 

D2I2C=(EB/R-2.*TM2(N) ) *EXP <- EB/ I R* TM? ( N) ) ) 

D2T2--D2T2A*D2T2B  +  D2T?C*D2T2D~?.*CT 

FUNIA=(  (SZKJ.I  1-SZ2I  J,  I  )  )*GA*SX1  (J, I ) *E A ) / ( R* < TM  I  (N ) **2 )  ) 

FUNIB»EXP(-EA/(R*TM1CN) ) ) 

FUN1C-ISZ2I J, I )*GB*SY1( J, I >*EB ) / ( R*( TM1 IN)**2)1 

FUN1C*EXPI-EB/IR*TMI(N) ) ) 

FUM1E=2.*CT*I2.*TM1 <N)-TM2(N)-T1M) 

FUN1»FUN1A*FUNIB+FUN1C*FUNID-F0N1E 

FUN2A=( (SZ3( J,I)-SZA( J, I ) )*GA*SX2(J, I ) *E A )/ I R*(TM2 ( N ! **2 ) ) 

FUM2a=EXPI-EA/(R*TM2{N) ) ) 

FUN2C=  I  S  Z>,  (J  ,  I  )  *SB*  SY2  IJ  .  I)  *EB  )  /  I  R*  (  TM?  (  N  !  *  *2  )  ) 

FUN2C=EXP(-E6/(R*TN2!N) ! ) 

FUN2E  =  2.*C"I*(TM1  (N1-TM2IN)  ) 

FUN2=FUN2A*FUN2B+FUN2C*FUN?D+FUN2E 

RHS1=TM1(N>*D1T1+TM2(N>*D1T?-FUN1 

RHS2»TM1(NJ*D2T1+TM2(N)*02T2-FUN2 

DErR=DlTlS02T2-DlT2*02Tl 

TMHN+1)* (RHS1*D2T2-RHS2*D1T2!/0ETR 

TW2(N+1)=(RHS2*DIT1-RHS1*D2T1)/DETR 


18  DELI=TK1(N+1)-TM1(N> 

19  DEL2=TM2(NU)-TM2<N) 
2C         TMNW1=TM1(N+1) 

21  TMNH2«TM2(N+11 

22  IF  (ABS(DELl).GT.O.l)  GO  TO  166 

23  IF  (ABS(OEL2).LF.0.1)  GO  TO  170 

24  166  COMTINUE 

25  168  FORf'AT(70X, 'TEMP  DID  NOT  CONVERGE  ',  2X ,  1  3  ,  •   OEl 1 = • , F9.3, 

l<   DEL2=',F9.3) 

26  PRINI    168,N,DELltDEL2 

27  170    TEHKJ,  I1=TMNW1 
2E  TEM2(J,l)=TMNW2 

25  IF     U.GE.3)    GO    TO    900 

3C  TMK  1)  =  TEM1  (  J,  I  ) 

31  TM2(ll=TEM2U.I) 

32  GO  TC  901 

33  9C0  TMlil  )---TEKlU-l,Hl> 

34  TM2(l)=TEM21J-li 1+1) 

35  901  CONTINUE 

36  PRINT  1C,,  I.ADVU  J.TEMKJ,  I),TEM2(J,I  1.PRF1  (  I  ) 
3?  165  CONTINUE 

C 

C      END  CE  ONE  ITERATION 

C 

38  300  CONTINUE 

C 

C      QUASILINPARI7.ATION  00  LOOP  ENDS  HERE 

C 
35  STOP 

4C  END 


41  SUBROUTINE  RKT  133 
C 

C      1HIS  SUBROUTINE  IS  USEO  TO  INTEGRATE  12  LINEARIZED 

C      EQUATIONS  S I  MULT ANEOSLY  BY  RUNGE  KUTTA  METHOD 

C 

42  COMMCN  X1,Y1,X2,Y2,AI,0,Z1,Z2,Z3,Z4,Z5,Z6,TP,TQ,DT,X0, 
1GA,GE,EA,R,AQ,V1,V2,Y0,EB,AN,CA,CI,AIM,CB,CC,CQ,C,AC, 
2CT,P,Wl,W2,W3,W4,TEr.l,TEM2,J,SAQ,SZ6 

43  DIMENSION    X 1  I  105 ) , Yl ( 105 ) , X2 ( 105 )  , Y2 ( 105 )  i A  I ( 1 05 ) , 

1Q( 105 ),Z1( 105 ),Z2 1105 ),Z3( 105), Z4( 105), Z5I105) ,Z6(105), 
2AK1C5),A2(105),A3(105),A4(105),B1(105),B2( 105 ) , B3 ( 105) , 
3B4dC5),CldC5),C2(105),C3(105),C4d05),Dld05),D2<105), 
4D3(lC5),D4d05),El(105),E2d05),E3(105),E4d05),Fl(105), 
5F2<1C5),F3(105),F4(105),C1(105),G2(105),G3(105),G4(105), 
6HKIC5) ,H2(105),H3d05),H41105),RK105),R2(105),R3(105), 
7R4(105),Sl(105),S2(l05),S3d05),S4(105),Tl(105),T2(105), 
8T3(1C5),T4I105),U1(105),U2(105),U3(105),U4(105),XO(105), 
9Y0UC0),TE(Ud0,102),TEM2d0,102J  tSAQ(10,102) ,SZ6I10,1U2) 
SP=AC/V1 
SQ=AC/V2 
DO    130    1=1,100 
TT1-R*TEM1 ( J-1,1 ) 
TT2=R*TEM2 ( J -  1 , I ) 
X0d)=Kl 
Y0(1 )=W2 
V=SAC(J-1,1> 
W=SZC(J-l, I ] 

Ald)=DT*(P*SP*XO(I)-SP*Xldl-GA*EXP(-EA/dl)*Xld)) 
Bld)  =  DT*(P*SP*YO(I)-SP*Yld)~GB*EXP(-EB/TTl)<Yld)  + 
1GA*EXP(-EA/TT1 ) *X  1  (  I  )  ) 
Cld)  =  DT*(SO*(Xl( I )-X2(I) )-GA*EXP(-EA/TT2)*X2( I) ) 
Did  )  =  DT*(SQ*(Y1( I)-Y2(I))-GB*EXP(-EB/Tr2)*Y21 I)+ 
1GA*EXP(-EA/TT2S*X2(I) ) 
Eld  )  =  DT*(AQ*Y2d  )-CQ*Qd)  ) 

FlA=tC*V-C*(V**2)/AN+V*K/(CA*AN)-W/(2.*CA)-W*(V**2)/ 
1(2.*CA*(AN**2) ) ) 
F2A=(C-2.*C*V/\N+W/(CA*AN)-W*V/(CA*(AN**2) ) ) 
F3A=(V/(CA*AN)-1./(2.*CA)-V**2/(2.*CA*(AN**2) ) ) 
rid )=DT*(P*F1A-P*V*F2A-P*W*F3A+Q! I)*F2A+Z6d)*F3A) 
Gld)  =  DT*(SP*ZU  I)-SQ*Z3(  I)  +  (ZH()-Z2(I)  ) *GA*EXP (-EA/TT 1 )  ) 
Hl(I)=DT*(SP*Z2d  )-SQ*Z4(I)  +  Z2(I  )  *EXP  (-EB/TT  1 )  *GB  ) 
Rid  )=DT*1S0*Z3(IH  (Z3d)-Z4d  )  > *GA*EXP (-EA/TT2 ) +P*AQ* ( CB-CC ) ) 
Sld)  =  DT*(SQ*Z4(  I  )  +  Z4d)*GB*EXP(-EB/TT2)-Z5d  )  *AQ-P*CC*AQ  1 
Tl(I)=DT*(P*2.*CI*AIM-2.*CI*AI (I ) ) 

U1A={AC*CQ-C*W+2.*C*V*W/AN+(W**2)*V/I2.*CA*(AN**2) ) 
1-W**2/(2.*CA*AN) ) 
U2A=(2.*C*W/AN+W**2/(2-*CA*(AN**2) ) ) 
U3A=(2.*C*V/AN-C+W*V/(CA*(AN**2) J-W/CCA*AN1 ) 
Uld )=DT*(CQ*Z5(I)+P*U1A-P*V*U2A-P*W*U3A+Q(I )*U2A+Z6! I) 
1*U3A) 

A2(I)=DT*(P*SP*X0( I )-SP*(Xl (IJ+A1 d )/2 . )-GA*EXP (-EA/TT1 ! 
1*1X11 I)+Al(I)/2.) ) 

B2(  I  )=Dr*(P*SP*Y0(  I  )-SP*(YK  I1+3U  I  )/2.  )-GB*EXP  ( -EB/TT  1 ) 
1*  CY1 CI ) +B1 d)/2. )+GA*EXP(-EA/TTl)*(Xld)+All I)/2. I) 
C2d)=0T*(SQ*<  (Xld)+AKI)/2.)-(X2d)+Cld)/2.)  !- 


1GA*EXP(-EA/TT2)*(X2(  l)+ClU>/2.  )  ) 

D2<I)=DT*<SQ*(<Yl(I)+Bl(I)/2.)-<Y2(I)+Dl(I)/2.))-GB*EX(> 
1I-EB/TT2)*(Y2(I )  +  Bl  <  I  ) /2.  ) +GA*EXP  (-EA/U2  J  *  i  X2  (  I  )  + 
2CKD/2.)) 

E2(I)=DT*IAQ*(Y2(I)+Dl(I)/2.)-CQ*(Q(I)+Flt I)/2.) ) 

F2I I )=DT*(P*FlA-P*V*F2A-Pn.'*F3A*(C( IJ+F1 (I )/2. )*F2A+ 
1(Z6< I)+U1( I)/2.)*F3A) 

G2(I)=DT*(SP*tZl( I )+GHI)/2. )-SQ*(Z3( Il+Rlll )/2.)+(( 
1ZH  I  )+Gl(  I  )/2.  )-(Z2(  I  )+Hl  (I  )/2.  )  )  *GA*EXP  (-EA/T1  1  )  ) 

H2(I)=DT*(SP*(Z2(I  >+Him/2.1-SQ*(Z4(I)  +  SllI)/2.  H 
1(Z2U)+H1( I)/2.)*EXP(-EB/TT1)*GB) 

R2(I  >=DT*(SQ*(Z3U  >+Rl  <  I  )/2.  )  + (  (  Z3  (  l) +R1  I  I  )/2.  )- (  Z'i  <  I  ) 
l+Sl(I)/2.) )*GA*EXP(-EA/TT2)+P*AQ*(CB-CC) ) 

S2(I )=DT*(SQ*(Z4( I)+Sl(I)/2.)+(Z4(I)+Sl( I)/2.)*GB* 
1EXP(-EB/TT2)-(Z5!I)4T1(I)/2.)*AQ-P*CC*AQ) 

T2(I  )=DT*(P*2.*CI*AII-!-2.*CI*IAI  (  I  )  +E1  (  I  )  /2.  )  ) 

U2( I)=DT*(C0*(Z5( I )+TKI)/2.)+P*UlA-P*V*U2A-P*W*U3A+ 
1(Q(  I  l+FKI  )/2.  )*U2A+[Z6<  I  )  +U1  (  I  )  12.  )*U3A  ) 

A3(I)=DT*(P*SP*X0(I J-SP*(X1(I)+A2( I ) /2 . )-GA*EXP I -EA/TT1 ) 
1*(X1( I)+A2(I)/2. J  J 

B3(I  )=DT*(P*SP*YO(  I  )-SP*(Yl< I)+B2(  I  )/2.  )-GB*EXP  (-EB/TT  I ) 
1*(Y1(I)+B2(I )/2. )+GA*EXP(-EA/TTl)*(Xl ( I ) * A2 ( I ) /2 . ) ) 

C3(1)=DT*(SQ*((X1(I )+A2(I)/2.)-(X2( I ) +C2 ( I ) /2. ) )- 
1GA*EXP(-EA/TT2)*U2U)+C2I  I)/2.)} 

D3(I)=DT*(SQ*((Y1(I  )+B2U)/2.)-(Y2U  ) +02  U)/2.  ))-GB*EXP 
1(-EB/TT2)*(Y2( I  )+B2U  1/2. > +GA*EXP (-EA/TT2 ) * < X2 ( I )+ 
2C2ID/2.)) 

E3(I)=DT*(AQ*(Y2(I ) 4D2 ( I )/2. )-CQ*( Q( I ) +F2 ( I  )/2. ) ) 

F3 ( I )=DT* ( P*F1A-P*V*F2A-P*W*F3A+ ( Q( I ) +F2 ( I ) /2. )*F2A+ 
1(Z6( I)+U2(I)/2.)*F3A) 

G3(I)=DT*(SP*(Z1(I )+G2(I)/2.  J-SQ*«Z3U)+R2(l)/2. )+( ( 
1Z1II )+G2(I )/2.)-(Z2(I)+H2(I)/2.)  ) *GA*EXP t-EA/TTl )  ) 

H3(  J  )  =  DT*(SP*(Z2(I>+!l2(I)/2.  )-SQ*(Z4C  I  )  +  S2l  I)  /2.  T+ 
1!Z2(  I1+H2I  I  )/2.)*EXP(-EB/TU)*GB) 

R3(I  )=t)T*(S0*(Z3(I)+!!2(I)/2.)  +  (  (  Z3  (  I  )  +  R2  <  [  1/2  .  )- (  Z<.  (  I  ) 
1  +  S2(  H/2.  )  )*GA*EXP(-EA/TT2)+P*AQ*(CB-CC)  ) 

S3 (I )  =  DT*(SQ*(ZMl )+S2( I)/2.)+(Z4( I )+S2( I)/2. )*GB* 
lEXP(-EB/rT2)-(Z5(I)+T2( I )/2. )*AQ-P*CC*AQ ) 

T3(I)  =  [)T*(P*2-*CI*AIM-2.*CI*(AI(I)  +  E2(I)/2.)) 

U3(I )=DT*(CQ*!Z5( I ) ♦ f2 ( I )/2. ) +P*U 1 A-P*V*U2A-P*W*U3A+ 
HQ(I)+F?H)/2.)*U2A+tZ6tI)+U2(I)/2.)*U3A) 

A'.(I)  =  DT*(P*SP*X0II)-S('*(X1(I)+A3(I))-GA*EXP(-EA/TT1)* 
HX1(I)  +  A3(I)  )) 

B'itI)  =  DT*(P*SP*YO(I  1-SP*1Y1(I)+B3(I  )  )-G8*EXP  (-EB/TT  1  )  * 
KYK  I  )+B3(I  )  )+GA*EXP(-EA/TTl)*(XKI  1+A31  I  )  )  ) 

C4(I)=DT*(SQ*I(X1(I)+A3(I) )- I X2 (I >+C3 ( I ) ) >-GA*EXP (-EA/ 
1TT2)*(X2(I )+C3(I  )  )  ) 

D'.(I)=DT*(SQ*(  (Y1II)  +  B3(I)  )-( Y2  (  I  ) +D3  ( I)  )  )-GB*EXP  ( -EB/ 
1TT2)*(Y2(  I  )+B3(I  )  )  +GA  -*EXP  (  -EA/TT2  )  *  (  X2  ( I  1+C3U  )  )  ) 

E4U)=DT*(AQ*(Y2i  11+03  1 [11-CQMQU )+F3(I) 1 1 

F4(I1=DT*(P*F1A-P*V*F2A-P*W*F3A+(Q(I)+F3  II ) )*F2A+(Z6(  I) 
1+U3I I) )*F3A) 

G4U)  =  DT*(SP*(Z1(I)+G3(1  )  )-S0*(Z3(I)+R3!  I)  )+(  (Zl(  IH 
lG3(I))-tZ2(I )+H3(I )>>*GA*EXP(-EA/TT1) ) 

H4(I)=DT*(SP*(Z2(I)+H3U) 1-SQ* ( Z4 ( I  1 +  S3 ( I > 1  +  ( Z2 ( I  1  + 
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C2 


i  i;  ;  i  i  ) 

RM  I) 

1))*GA 

S4(I) 

ld2)- 

T4(  I  > 

114(1) 

HF3I  r 

XIII* 

Yl(  U 

X2(  [  + 

Y2(  1  + 

aiih 

0(1  +  1 
2111  + 
Z2(d 
Z3II  + 
Z4  (I  4 
1 5  (  1  4 
Z61I4 
130  CONTI 
RETUR 
END 


>*EXP 
=  DT*  t 

*EXP( 

=or*( 

(Z5(I 
=  DT*( 
=  DT*( 
)  >*U2 
1)  =  X1 
l)  =  Yl 
1)=X2 
1)  =  Y2 

I  )  =  AI 
)  =  0d 
1)  =  Z1 

II  =  22 
1)  =  Z3 
l)  =  Z4 
1)  =  Z5 
1)  =  Z6 
NUE 

N 


(-E8 

SQ*( 
-EA/ 
SQ*( 
I+T3 
P*2. 
C0*( 
A4(Z 

(I)  + 
(I  )4 
(1)4 
(I  )4 
(1)4 
)4l. 
(1)4 
(1  )4 
(1)4 
(I  )  + 
(I  )4 
(  I  )4 


I  1  ) 
(I) 
2)4 
(  I  ) 
)  )* 

I*A 
(  I  ) 
I  H 
fb. 
fb. 
fb. 
fb. 
fb. 
*( 
fb. 
fb. 
fb. 
/b. 
fb. 
fb. 


B) 

3d)  ) 
AQ*(C 
3(1)  ) 
-P*CC 
2.*C 
3(  I)  ) 
(I)  )* 
All  I) 
bid) 
C1U) 
Did) 
Eld) 
(1)42 
GUI  ) 
Hid) 
Rl  (I) 
SKI) 
Till) 
Ul  d  ) 
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4( (Z3I I 
B-CC)  ) 
4(Z'.d) 
*AQ) 
I*(AI(I 
4P*U1A- 
U3A) 
42.*A2( 
42.*B2( 
42.*C2( 
42.*D2( 
42.*E2( 
*F2(  I) 
42.*G2( 
42.*H2( 
42.*R2( 
42.*S2( 
42.*T2( 
42.*U2( 


)4R3( I  )  )-(ZM  I )4S3(I ) 
4S3( I ) )*GB*EXP(-EB/ 


)4E3d))) 
P*V*U2A-P*W*U3A+tQl  I ! 


1)42. *A3( 
1  )42.*B3( 
I  )42.*C3( 
1)42. *D3( 
I  )42.*E3( 
42.*F3(I) 
1  )42.*G3( 
I)+2.*H3( 
I  )42.*R3( 
1)42. *S3( 
I  I  42. *T 3  1 
1 )42.*U3( 


I  )+M{  I 
I  )4B4( I 
1  )4CAd 
I  )4DMI 
I)4E<.d 
4FMII  ) 
I)4GMI 
1 )4H4(  I 
I )4R4( I 
I)4S4d 
I  )4T4d 
IdUM  I 
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C 

C      THIS  SUBROUTINE  IS  USED  TO  INVERT  A  SIX  BY  SIX 
C      MATRIX  ENCOUNTERED  IN  THE  CALCULATION  OF  SIX 
C      INTEGRATION  CONSTANTS.  THIS  IS  SUPPLIED  BY  IBM. 
C 

DIMENSION  A(1),B(1! 
C 

C         FORWARD  SOLUTION 
C 

TOL=C.O 

KS=0 

JJ=-N 

DO  65  J=1,N 

JY=Jtl 

JJ  =  JJ  +  IJ+1 

BIGAsO 

IT=JJ-J 

DO  30  I=J,N 
C 

C         SEARCH  FOR  MAXIMUM  COEFFICIENT  IN  COLUMN 
C 

IJ= IT + I 

IFtACS(BICA)-ABS(A( IJ) ) )  20,30,30 
20  BIGA=A(IJ) 

IMAX=I 
30  CONTINUE 
C 

C         TEST  FOR  PIVOT  LESS  THAN  TOLERANCE  (SINGULAR  MATRIX) 
C 

IF(APS(BIGA)-TOL)  35,35,'i0 
3b  KS=1 

RETURN 
C 

C         INTERCHANGE  ROWS  IF  NECESSARY 
C 

40  Ii=J*N*!J-2) 

IT=I^AX-J 

DO  5C  K=J,N 

1 1=1  UN 

12=1  HIT 

SAVE=A( II) 

A(  I1)=A(I2) 

A(I2)=SAVE 
C 

C         DIVIDE  EQUATION  BY  LEADING  COEFFICIENT 
C 

50  A(II)=A(I1)'BIGA 

SAVE=BI (MAX) 

B(IMAX)=B(J) 

Bl J)=SAVF/BIGA 
C 

C         ELIMINATE  NEXT  VARIABLE 
C 

IFIJ-N)  55,70,55 


55  IQS=N*(J-1) 

DO  65  IX=JY,N 

IXJ=IQS+IX 

IT=J-IX 

DO  6C  JX=JY,N 

IXJX=N*UX-1)  +  IX 

JJX=IXJX+IT 
60  A(IXJX)=M IXJX)-(A(IXJ)*A(JJX) ) 
65  B(1X)  =  B(  IX)-(B(J)*AUXJ)  > 
C 

C         BACK  SOLUTION 
C 

70  NY=N-1 

IT=N*N 

DO  8C  J=1,NY 

ia=ii-j 

1B=N-J 
IC=N 

DO  8C  K=1,J 

BUB)=B(IB)-A(IA)*BIIC) 
IA«IA-N 
80  IC=IC-1 
RETURN 
END 
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The  importance  of  quantitative  techniques  in  decision  making  empha- 
sizes the  need  of  efficient  techniques  as  a  tool  for  solving  management 
problems.   However,  fairly  powerful  algorithms  are  not  yet  available  in 
solving  dynamic  management  problems  involving  differential  equations. 

The  two  point  boundary  value  problem  with  non-linear  differential 
equations  provides  such  an  example.   The  nonlinearity  in  the  performance 
equations  does  not  allow  the  application  of  superposition  principle. 

Quasilinearization  helps  overcome  this  difficulty.   It  linearizes 
the  non-linear  equations  and  provides  an  algorithm  which  would  give  the 
solution  by  an  iterative  procedure. 

The  purpose  of  this  work  is  to  investigate  the  effectiveness  of  this 
recently  developed  tool  in  solving  various  industrial  management  problems. 

First  a  brief  introduction  and  computational  procedure  of  quasilinear- 
ization is  given.   Then  its  application  to  an  advertisement  problem  with 
two  state  variables  and  one  control  variable  is  discussed  in  detail. 

Next  is  discussed  the  application  of  quasilinearization  to  an  adver- 
tisement and  production  problem.   This  model  has  six  state  variables  and 
three  control  variables.   In  addition,  the  profiles  are  fairly  unstable 
due  to  the  rapid  change  of  variables  with  time. 

It  was  concluded: 

1.  Choosing  the  initial  approximations  to  start  the  solution  is  not 
difficult  in  most  cases. 

2.  The  convergence  rate  is  almost  independent  of  the  choice  of  the 
initial  approximations. 

3.  This  algorithm  converges  quadratically ,  if  it  does  converge. 


4.  For  rapidly  increasing  or  rapidly  decreasing  profiles,  first 
variational  and  second  variational  techniques  seemed  to  have 
failed.   On  the  other  hand,  this  method  encountered  no  problem 
in  converging  to  the  optimal  solution. 

5.  Because  of  the  intimate  association  between  the  boundary  value 
problems  and  the  optimization  and  control  problems,  this  technique 
may  provide  a  useful  tool  for  the  systems  analysts. 


